REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collertion  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  AND  SUBTITLE  /  /  /  /)  /) 


3.  REPORT  TYPE  AND  DATES  COVERED 


6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADORESS(ES) 
DEPARTMENT  OF  THE  AIR  FORCE 
AFTT/CI 

2950  P  STEET,  BLDG  125 
WRIGHT-PATTERSON  AFB  OH  45433-7765 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  118.  SECURITY  CLASSIFICATION  119.  SECURITY  CLASSIFICATION  I  20.  LIMITATION  OF  ABSTRACT 

OF  REPORT  I  OF  THIS  PAGE  I  OF  ABSTRACT  I 


NSN  7540-01-280-5500 


DXxo 


XjLiliX  1 


Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 
298-102 


GENERAL  INSTRUCTIONS  FOR  COMPLETING  SF  298 

The  Report  Documentation  Page  (RDP)  is  used  in  announcing  and  cataloging  reports.  It  is  important 
that  this  information  be  consistent  with  the  rest  of  the  report,  particularly  the  cover  and  title  page. 
Instructions  for  filling  in  each  block  of  the  form  follow.  It  is  important  to  stay  within  the  iines  to  meet 
optical  scanning  requirements. 


Block  1.  Agency  Use  OnW  (Leave  blank). 


Block  2.  Report  Date.  Full  publication  date 
including  day,  month,  and  year,  if  available  (e.g.  1 
Jan  88).  Must  cite  at  least  the  year. 

Blocks.  Type  of  Report  and  Dates  Covered. 


State  whether  report  is  interim,  final,  etc.  If 
applicable,  enter  inclusive  report  dates  (e.g.  10 
Jun87-30  Jun88). 

Block  4.  Title  and  Subtitle.  A  title  is  taken  from 


the  part  of  the  report  that  provides  the  most 
meaningful  and  complete  information.  When  a 
report  is  prepared  in  more  than  one  volume, 
repeat  the  primary. title,  add  volume  number,  and 
include  subtitle  for  the  specific  volume.  On 
classified  documents  enter  the  title  classification 
in  parentheses. 

Blocks.  Funding  Numbers.  To  include  contract 


and  grant  numbers;  may  include  program 
element  number(s),  project  number(s),  task 
number(s),  and  work  unit  number(s).  Use  the 
following  labels; 


Contract 

PR  - 

Project 

Grant 

TA  - 

Task 

Program 

WU  - 

Work  Unit 

Element 

Accession  No 

Blocks.  Author(s).  Name(s) of person(s) 
responsible  for  writing  the  report,  performing 
the  research,  or  credited  with  the  content  of  the 
report.  If  editor  or  compiler,  this  should  follow 
the  name(s). 

Block?.  Performing  Organization  Name(s)  and 


Address(es).  Self-explanatory. 

Block  8.  Performing  Organization  Report 


Number.  Enter  the  unique  alphanumeric  report 
number(s)  assigned  by  the  organization 
performing  the  report. 

Block  9.  Sponsor! no/Monitori no  Agency  Name(s 


and  Address(es).  Self-explanatory. 

Block  10.  Sponsoring/Monitoring  Agency 
Report  Number.  (If  known) 

Block  11.  Supplementary  Notes.  Enter 
information  not  included  elsewhere  such  as: 
Prepared  in  cooperation  with...;  Trans,  of...;  To  be 
published  in....  When  a  report  is  revised,  include 
a  statement  whether  the  new  report  supersedes 
or  supplements  the  older  report. 


Block  12a.  Distribution/Availabilitv  Statement. 
Denotes  public  availability  or  limitations.  Cite  any 
availability  to  the  public.  Enter  additional 
limitations  or  special  markings  in  all  capitals  (e.g. 
NOFORN,  REL,  ITAR). 


See  DoDD  5230.24,  "Distribution 
Statements  on  Technical 
Documents." 

See  authorities. 

See  Handbook  NHB  2200.2. 
Leave  blank. 


DOE  - 
NASA- 
NTIS  - 


Block  12b.  Distribution  Code. 


NASA 

NTIS 


Leave  blank. 

Enter  DOE  distribution  categories 
from  the  Standard  Distribution  for 
Unclassified  Scientific  and  Technical 
Reports. 

Leave  blank. 

Leave  blank. 


Block  13.  Abstract.  Include  a  brief  (Max/mum 
200  words)  factual  summary  of  the  most 
significant  information  contained  in  the  report. 

Block  14.  Subject  Terms.  Keywords  or  phrases 
identifying  major  subjects  in  the  report. 

BlockIS.  Number  of  Pages.  Enterthe  total 


number  of  pages. 

Block  16.  Price  Code.  Enter  appropriate  price 
code  (NTIS  only). 

Blocks  17.  - 19.  Security  Classifications.  Self- 


explanatory.  Enter  U.S.  Security  Classification  in 
accordance  with  U.S.  Security  Regulations  (i.e., 
UNCLASSIFIED).  If  form  contains  classified 
information,  stamp  classification  on  the  top  and 
bottom  of  the  page. 

Block  20.  Limitation  of  Abstract.  This  block  must 


be  completed  to  assign  a  limitation  to  the 
abstract.  Enter  either  UL  (unlimited)  or  SAR  (same 
as  report).  An  entry  in  this  block  is  necessary  if 
the  abstract  is  to  be  limited.  If  blank,  the  abstract 
is  assumed  to  be  unlimited. 


★  U.S.GPO:  1 993-0-336-043 


Standard  Form  298  Back  (Rev.  2-89) 


A  Finite  Element  Model  for  Harmonic  Response  of  a 
Viscoelastic  Sandwich  Beam 


A  Thesis 


Presented  to 
The  Faculty  of  the 

School  of  Engineering  and  Applied  Science 
University  of  Virginia 


In  Partial  Fulfillment 
of  the  Requirements  for  the  Degree 
Master  of  Science  in  Civil  Engineering 

by 

Richard  Allen  Maddox 
May  1996 


19961212  021 


HI  ClAIMffi  NOm 


TfflS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE  COPY 
FURNISHED  TO  DUC  CONTAINED 
A  SIGNHTCANT  NUMBER  OF 
COLOR  PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY  ON  BLACK 
AND  WHITE  MICROnCHE. 


APPROVAL  SHEET 


The  thesis  is  submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


Master  of  Science  in  Civil  Engineering 


This  thesis  has  been  read  and  approved  by  the  Examining  Committee: 


Thomas  T.  Baber,  Thesis  Advisor 


,  Faculty  Committee 


Accepted  for  the  School  of  Engineering  and  Applied  Science; 


Dean,  School  of  Engineering  and 
Applied  Science 


May  1996 


1 


TABLE  OF  CONTENTS 


Page 

LIST  OF  FIGURES  iii 

LIST  OF  TABLES  iv 

LIST  OF  SYMBOLS  v 

ABSTRACT  vii 

1.  INTRODUCTION  I 

1.1  Industrial  Vibration  Absorption  1 

1 .2  Sandwich  Beams  for  TMD  Applications  4 

1.3  Viscoelasticity  5 

2.  SANDWICH  BEAM  MODELS  8 

2.1  Single  Variable  Models  8 

2.2  Multiple  Variable  Models  14 

3.  FINITE  ELEMENT  FORMULATION  FOR 

DAMPED  SANDWICH  BEAMS  24 

3.1  Prior  Research  24 

3.2  Assumptions  26 

3.3  Assuming  Shape  Functions  33 

3.4  FORTRAN  Program  38 

3.5  Equation  of  Motion  38 

3.6  FORTRAN  Solution  Process  39 

4.  NUMERICAL  STUDIES  41 

4.1  Test  Cases  41 

5.  CONCLUSIONS  49 

BIBLIOGRAPHY  51 

APPENDICES 

1.  Mass  and  Stiffness  Matrices  A1 

2.  FORTRAN  Program  A2 

3.  Ansys  Program  A3 


ii 


LIST  OF  FIGURES 


#  Page 

1  Typical  Tuned  Mass  Damper  2 

2  Beam-type  TMD  3 

3  Calculated  depth,  d,  compared  with  Actual  depth  12 

4  Bai  and  Sun  Shear  Deformation,  T  ,15 

5  Sandwich  Beam  Construction  1 5 

6  Sandwich  Beam  Displacement  Field  19 

7  Reducing  Displacement  Field  20 

8  Plot  showing  effect  of  core  thickness  on  shear  deformation,  32 

9  1st  Order  Shape  Functions  35 

10  2nd  Order  Shape  Functions  36 

11  Displacement  Vector  37 

12  Specimen  1  Test  Data  44 

13  Lu  and  Douglas  Test  Data  45 

14  Specimen  2  Test  Data  46 

1 5  Ansys  Model  48 


iii 


LIST  OF  TABLES 


# 

1  Test  Case  Characteristics 


2  Ansys  Results 


LIST  OF  SYMBOLS 


Oy  matrix  element  of  the  ith  row,yth  column 
A  cross  sectional  area 

c  viscous  damping 

d  depth  of  shear  block 

e  2(1 +v,) 

E  Young's  modulus 

complex  Young's  modulus  of  core 
Ex  Young's  modulus  for  top  layer 

£2  Young's  modulus  for  bottom  layer 

E*  complex  Young's  modulus 

/  natural  frequency 

F  load  vector 

F(j  harmonic  loading  vector 

G  shear  modulus 

complex  shear  modulus  of  core 
G*  complex  shear  modulus 

h  core  layer  thickness 

/ii  top  layer  thickness 

h2  bottom  layer  thickness 

i  square  root  of  - 1 ,  element  node 

/  moment  of  inertia 

Imp  driving  point  mechanical  impedance 
j  square  root  of  - 1 ,  element  node 

k  stiffness,  element  node 

k*  complex  stiffness 

K  stiffness  matrix 

K*  complex  stiffness  matrix 

L  beam  length 

m  mass 

M  cross  sectional  moment, 

M  mass  matrix 

N  shape  function 

p  transverse  load 

q  external  load 

Q  harmonic  loading 

S  sum  of  shear  forces 

t  time 

T  Kinetic  energy 

u  longitudinal  displacement 


V 


w*  longitudinal  displacement  of  bottom  face 

Wo  longitudinal  displacement  at  centroid  of  core 

u  longitudinal  displacement  of  top  face 

w,  longitudinal  displacement  of  top  layer  at  centroid 
W2  longitudinal  displacement  of  bottom  layer  at  centroid 

w  longitudinal  acceleration 

U  displacement  vector,  potential  energy 

U  acceleration  vector 

w  transverse  displacement 

■w’  transverse  displacement  of  bottom  face 

-w^  transverse  displacement  of  core 

Wq  transverse  displacement  at  centroid  of  core 

w  transverse  displacement  of  top  face 

W]  transverse  displacement  of  top  layer  at  centroid 

W2  transverse  displacement  of  bottom  layer  at  centroid 

>V|'  rotation  displacement  of  top  face 

W2  rotational  displacement  of  bottom  face 

X  Cartesian  coordinate 

y  Cartesian  coordinate 

z  Cartesian  coordinate 

a  X  dependent  portion  of  longitudinal  displacement,  w 

p  X  dependent  portion  of  transverse  displacement,  w 

6  damping  factor 

y  shear  strain 

O  transverse  normal  deformation  in  core 

V  Poisson's  ratio 

Vc  Poisson's  ratio  of  core 

p  mass/unit  length 

Pi  mass/unit  length  of  top  layer 

P2  mass/unit  length  of  bottom  layer 

T  shear  force 

CO  forcing  frequency 

Q  forcing  frequency 

4^  shear  deformation 


VI 


ABSTRACT 


After  a  discussion  of  the  possible  use  of  the  sandwich  beam  design  as  a  means  of 
incorporating  damping  into  tuned  mass  dampers,  the  development  of  a  finite  element 
model  for  viscoelastic  sandwich  beam  analysis  is  presented. 

The  requirement  for  future  use  of  the  viscoelastic  sandwich  beam  design  in  TMDs  is 
to  have  a  method  by  which  the  response  of  a  beam  with  any  given  dimensions  or  end 
conditions  can  be  determined.  A  survey  of  previous  published  research  into  sandwich 
beam  behavior  demonstrates  the  need  for  the  finite  element  approach.  The  theory  and 
assumptions  presented  in  this  previous  research,  provide  the  basis  for  the  development  of 
the  approximation  model  presented  here.  A  finite  element  model  is  presented  that  utilizes 
sandwich  beam  theory  developed  for  thick  damping  layers.  The  model  is  constructed 
using  standard  beam  and  bar  shape  ftmctions.  Through  an  approximation  introduced  by 
observing  the  spatial  variation  of  the  core  shear  deformation,  it  was  possible  to  eliminate 
all  core  variables  and  express  the  element  behavior  in  terms  of  nodal  displacements  of  the 
top  and  bottom  face  plates  only.  A  FORTRAN  computer  program  was  written  to  perform 
the  finite  element  analysis. 

Using  experimental  data  from  previous  research  and  data  obtained  through  use  of 
computer  analysis  software,  results  from  this  model  are  compared.  The  results  are  used 
to  support  the  conclusion  that  this  working  model  can  be  used  in  the  future  for  the 
analysis  of  viscoelastic  sandwich  beams. 


vii 
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1.0  INTRODUCTION 

1.1  Industrial  Vibration  Absorption 

When  industrial  machinery  is  located  within  a  structure,  the  resulting  dynamic  forces 
can  have  severe  effects  on  the  structure.  The  repetitive  forces  of  many  industrial 
machines  may  induce  excessive  floor  vibrations.  Over  time  the  repeated  back  and  forth 
deflection  of  floor  members  can  cause  fatigue  and  potential  failure.  If  the  frequency  of 
the  machinery  coincides  with  the  structure's  natural  frequency,  the  effects  are  most 
severe.  In  addition  to  possible  damage  to  the  structure,  the  machinery  producing  the 
disturbance  moves  with  the  floor  to  which  it  is  attached.  Therefore  damage  to  or 
malfunction  of  the  machine  itself  becomes  a  concern,  especially  in  modem  machinery 
containing  sensitive  computer  components.  An  example  of  this  occurred  when  the 
upgrade  of  machinery  in  one  factory  led  to  severe  difficulties.  The  faster  operating 
speeds  of  the  new  equipment  induced  a  large  vibration  in  the  floor.  In  response  to  the 
motion  of  the  floor,  the  machinery  automatically  shut  itself  down  as  a  result  of  loss  of 
alignment. 

To  improve  a  structure’s  behavior  under  periodic  loading  there  are  several  alternatives. 
The  structure  can  be  stiffened  to  raise  the  resonant  frequency,  mass  can  be  added  to  the 
stmcture  to  lower  the  resonant  frequency,  isolating  pads  can  be  used  to  reduce  the  transfer 
of  energy  from  the  machinery  to  the  structure,  or  an  energy  absorbing  device  may  be 
installed.  Stiffening  the  structure  or  increasing  its  mass  will  shift  its  natural  frequency 
out  of  the  operating  range  of  the  machinery,  but  may  require  considerable  modifications 
to  the  structure.  This  would  mean  increased  costs  to  new  consfruction  and  expensive 
alterations  to  an  existing  stmcture  which  may  not  be  physically  possible  (Auerbach, 

1995).  The  use  of  isolating  pads  can  be  beneficial,  but  will  not  always  alleviate 
vibrations  problems  alone,  especially  when  the  vibrations  are  large. 

For  the  purpose  of  absorbing  vibrations,  active  and  passive  techniques  can  be 
employed  which  would  dissipate  the  energy  being  applied  to  the  stmcture  thus  reducing 
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the  induced  vibration.  Passive  devices  tend  to  be  more  economical  and  reliable  than 
active  devices  which  depend  upon  intricate  electronics.  One  technique  for  vibration 
control  is  that  of  installing  tuned  mass  dampers  (TMDs)  on  the  structure.  A  TMD  is  a 
passive  vibration  absorber  which  dissipates  a  portion  of  the  energy  from  a  structure  to 
which  it  is  attached  by  vibrating  independently  of  the  structure.  The  natural  frequency  of 
the  TMD  is  "tuned"  (by  design)  to  an  optimum  natural  frequency  allowing  it  to  absorb  the 
greatest  amount  of  energy  when  the  structure  to  which  it  is  attached  reaches  resonance. 


Figure  1  -  Typical  Tuned  Mass  Damper 

TMDs  are  effectively  used  in  mechanical  engineering  systems,  including  machinery, 
automotive  and  aircraft  engines,  and  ships,  and  to  reduce  wind  induced  vibration  of  tall 
buildings  and  structures  (Xu,  Kwok,  and  Samali,  1992).  The  design  of  TMDs  involves 
the  consideration  of  its  three  components:  mass,  stiffness,  and  damping.  The  greater  the 
mass  of  the  TMD,  the  greater  the  amount  of  energy  it  will  absorb.  The  amount  of  mass 
which  can  be  used  however  is  limited  by  practical  concerns.  Adding  too  much  weight 
would  require  modification  of  the  structure  to  support  it.  Size,  cost,  and  constructability 
all  dictate  that  the  TMD’s  mass  be  small  relative  to  the  floor  to  which  it  is  attached 
(Auerbach,  1995). 
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Since  only  small  variations  in  the  mass  are  possible,  the  stiffness  of  the  TMD 
(represented  as  k  in  Figure  1)  becomes  critical  in  determining  the  natural  frequency  of  the 
TMD.  The  required  stiffness  may  be  obtained  in  the  design  of  the  TMD  through 
geometry,  material  composition,  and  cormection  methods.  In  general,  increasing  the 
stiffness  of  the  TMD  will  decrease  the  amount  of  damping  (represented  as  c  in  Figure  1) 
present  and  vice  versa  (Bai  and  Sun,  1995).  The  design  of  TMDs  requires  balancing 
these  components  to  absorb  the  greatest  amount  of  vibratory  energy. 

Specifically  the  device  which  this  research  is  concerned  with  is  the  beam-type  TMD. 
The  use  of  a  beam-type  design  allows  for  construction  of  a  device  which  can  attach  easily 
to  the  floor  system  in  a  structure  and  would  be  compact  enough  to  fit  under  the  floor 
without  major  modifications  to  the  structure.  A  beam-type  TMD  would  be  equipped  with 
movable  masses  at  each  end  for  fine  tuning  the  device  during  installation. 


Mass  Mass  Mass 


Figure  2  -  Beam-type  TMD 


The  disadvantage  of  the  beam  design  is  that  the  one  mode  of  vibration  possible  in  the 
lumped  mass  design  in  the  Figure  1  is  now  increased  to  many  possible  modes.  Carefiil 


analysis  must  accompany  any  design  to  ensure  that  the  proper  mode  of  the  beam-type 
TMD  is  excited. 
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Previous  research  at  the  University  of  Virginia  has  produced  an  optimal  design 
method  for  beam-type  TMDs.  To  test  the  accuracy  of  his  design  method,  a  beam-type 
TMD  was  constructed  in  the  lab  by  Auerbach  (1995).  A  model  floor  constructed  by 
Sgambati  (1994)  was  used  as  the  platform  for  the  experiment.  Walsh  (1994)  added 
braces  to  the  floor  and  determined  the  natural  frequency  and  damping  of  the  braced  floor. 
Using  this  data,  Auerbach  was  able  to  design  a  beam-type  TMD  for  the  braced  floor.  He 
constructed  such  a  TMD  and  attached  it  under  the  floor  system.  As  predicted,  the  TMD 
had  the  stiffness  necessary  to  affect  the  floor's  response  to  applied  loads.  What  was 
missing  however  was  the  damping  component  of  the  TMD.  Commercially  made 
pneumatic  dashpots  and  attaching  butyl  rubber  strips  to  the  outer  faces  of  the  TMD  beam 
were  able  to  provide  a  small  amount  of  damping  to  the  system.  However,  the  optimal 
damping  needed  could  not  be  reached. 

1.2  Sandwich  Beams  for  TMD  Applications 

The  problem  then  becomes  to  develop  a  beam-type  TMD  which  has  stifi&iess  to  resist 
deflection  and  yet  still  provides  damping  of  the  vibration.  For  that  purpose,  this  research 
was  vmdertaken  into  the  use  of  the  sandwich  beam  design.  Plantema  (1966)  defines 
sandwich  construction  as  a  layer  type  construction,  consisting  of  thin  sheets  of  high- 
strength  material  between  which  a  thick  layer  of  low  average  strength  and  density  is 
sandwiched.  Although  multiple  layers  are  possible,  this  research  is  only  concerned  with 
the  simple  three-layer  sandwich  beam.  The  two  thin  sheets  are  called  the  faces,  and  the 
intermediate  layer  is  the  core  of  the  sandwich. 

In  the  aerospace  industry,  where  the  use  of  sandwich  construction  is  most  prevalent 
mainly  due  to  its  weight  savings  and  ease  of  mass  production,  the  material  of  the  faces 
may  be  aluminum  alloy,  reinforced  plastic,  titanium,  heat-resistant  steel,  etc.  For  the 
core  the  material  and  the  geometric  shape  can  vary  widely.  Some  of  the  main  types  of 
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cores  used  today  are  honeycomb  cores,  corrugated  cores,  balsa  wood,  and  expanded 
materials  such  as  cellulose  acetate  or  synthetic  rubber  (Plantema,  1966). 

1.3  Viscoelasticity 

To  achieve  a  high  level  of  damping  from  the  TMD,  a  core  of  viscoelastic  material 
would  be  used.  Viscoelasticity  is  the  property  of  special  materials  such  as  rubber,  which 
exhibit  characteristics  of  being  viscous,  with  stresses  proportional  to  the  rate  of  strain,  and 
of  being  elastic,  with  stresses  proportional  to  the  strains.  More  generally,  viscoelastic 
materials  are  those  materials  whose  constitutive  law  relates  the  stresses  and  their  time 
derivatives  to  the  strains  and  their  time  derivatives. 

When  analyzing  isotropic  elastic  materials,  the  relationships  between  stress  and  strain 
that  a  material  exhibits  can  be  represented  by  three  values:  the  Young's  modulus,  E,  the 
Shear  modulus,  G,  and  the  Poisson's  ratio,  v.  For  isotropic  elastic  materials  where 
stresses  are  applied  in-plane  with  the  material  (plane  stress),  the  Young's  modulus  is 
related  to  the  shear  modulus  by  the  following  relationship: 


E 

2(1 +  v) 


[1] 


The  behavior  of  viscoelastic  materials  is  not  as  simple  to  define  as  that  of  elastic 
material.  As  stress  is  applied  to  a  viscoelastic  material,  it  begins  to  deform  elastically, 
but  as  the  stress  continues  over  time,  there  is  a  certain  amount  of  mechanical  loss  in  the 
material  caused  by  the  strain.  For  linearly  viscoelastic  materials,  a  general  relationship 
between  uniaxial  stress  and  strain  is  given  by  the  differential  equation: 


«o  +«i 


dt 


'  +  Cl^ 


dV 


■+a 


"  dt") 


(u  u  ^  u 

I  °  'a/  ' ar'  at" 


[2] 
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For  general  loading  it  is  necessary  to  specify  the  response  in  terms  of  this  differential 
equation,  leading  to  a  transient  response. 

In  the  case  of  response  of  a  structure  to  harmonic  excitation,  a  somewhat  simpler 
approach  is  possible.  One  early  researcher  to  use  this  approach  was  E.M.  Kerwin  Jr. 
(DiTaranto,  1965).  If  the  excitation  is  taken  in  complex  form  as  ,  then  the  stresses 
and  strains  may  both  be  assumed  to  be  of  the  same  form.  That  is 


a  =  [3] 

8  = 

Substituting  these  equations  into  the  differential  constituitive  equation  [2],  one  obtains 
[a^+JCla^  +  a2+--HJ^y  a„]G  =  [&„  +  yQb,  +  (70)^62 +--+(7Q)^6p]8 

[4] 

This  suggests  that  the  stress-strain  relationship  under  harmonic  forcing  may  be  expressed 
in  the  complex  form 

a  =  E(l  +  j5,)e  [5] 

The  analysis  of  viscoelastic  materials  under  harmonic  forcing  therefore  requires  the 
use  of  complex  moduli  to  represent  both  the  elastic  response  of  the  material  and  the 
response  due  to  the  viscous  dissipation  of  energy  over  time.  The  fact  that  the  ratio  of 
stress  to  strain  is  a  complex  quantity  signifies  that  strain  lags  in  phase  behind  stress  by  an 
angle  the  tangent  of  which  is  called  the  loss  or  damping  factor.  The  resulting  strain  of 
viscoelastic  material  under  a  sufficiently  small  applied  stress  is  a  function  of  time  alone 
and  not  of  stress  magnitude  (Snowdon,  1968). 


In  many  problems  of  interest  much  of  the  energy  dissipation  occurs  in  shear, 
following  the  above  reasoning,  Snowdon  describes  a  complex  shear  modulus,  G*,  to  be 
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used  in  the  analysis  of  viscoelastic  materials: 

G*  =  complex  shear  modulus  =  G  (1  +76)  [6] 

where  j  =  sq  root  of  -1 

G  =  the  dynamic  shear  modulus  =  is  the  real  part  of  the  modulus 
5  =  the  loss  or  damping  factor  =  is  the  ratio  of  the  real  part  to  the  imaginary  part 
G6  =  the  imaginary  part  =  is  a  measure  of  the  mechanical  loss  associated  with  the 

shear  deformation  of  the  material 

G*  and  G  are  functions  of  temperature  and  of  the  angular  frequency  at  which  the  stress 
and  strain  vary  sinusoidally  with  time.  As  temperature  decreases  and  as  frequency 
increases  the  shear  modulus  increases.  The  increase  continues  towards  a  point  where  the 
dynamic  modulus  becomes  so  large  that  the  characteristic  resilience  of  a  rubberlike 
material  is  no  longer  apparent.  An  inextensible  (or  glasslike)  state  is  reached  (Snowdon, 
1968).  This  point  is  known  as  transition  and  the  frequency  and  temperature  at  which  this 
occurs  is  the  transition  frequency  and  the  transition  temperature.  Viscoelastic  material 
with  a  low  transition  frequency  will  have  a  maximum  damping  factor  occurring  at  a  low 
frequency.  Those  materials  which  produce  high  damping  will  have  low  transition 
frequencies.  Examples  of  high  damping  viscoelastic  materials  are  plasticized  polyvinyl 
butyral  resin,  Thiokol  RD,  plasticized  polyvinyl  acetate,  and  butyl  rubber  (filled  with  40 
parts  by  weight  of  medium  processing  channel  (MFC)  black  per  100  parts  of  rubber) 
(Snowdon,  1968). 

Similar  to  the  complex  shear  modulus  is  the  complex  Young's  modulus: 

E*  -  complex  Young's  modulus  =  E  (1  +  jh)  [7] 

The  Poisson's  ratio  for  viscoelastic  materials  approaches  the  perfectly  elastic  state  for 
which  the  ratio  is  1/2.  For  calculations  in  this  paper,  a  Poisson's  ratio  of  0.49  will  be 
used.  For  purposes  of  this  research,  the  complex  Yoimg’s  modulus  is  assumed  to  be 
related  to  the  complex  shear  modulus  by  equation  [1]. 
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2.0  SANDWICH  BEAM  MODELS 

2.1  Single  Variable  Models 

Although  the  use  of  the  viscoelastic  sandwich  beams  have  not  been  previously 
researched  specifically  for  the  use  in  tuned  mass  dampers,  they  have  been  the  object  of 
investigation  in  the  reduction  of  the  transmission  of  vibratory  energy  in  structures  and 
mechanical  systems,  especially  in  the  aerospace  industry,  for  many  years. 

One  of  the  first  attempts  to  formulate  static  and  dynamic  bending  equations  for 
viscoelastic  sandwich  beams  was  done  in  a  paper  by  R.A.  DiTaranto  (1965).  Using  some 
basic  assvunptions  set  forth  by  Kerwin,  DiTaranto  was  able  to  develop  a  sixth-order, 
complex,  homogeneous  differential  equation  for  longitudinal  displacements.  He  was 
then  able  to  show  how  the  solution  of  this  equation,  once  boundary  conditions  were 
satisfied,  would  yield  the  natural  frequencies  and  damping  factors  for  viscoelastic 
sandwich  beams. 

The  assumptions  which  DiTaranto  followed  were; 

1.  For  the  beam  cross  section,  there  is  a  neutral  axis,  whose  location  varies  with 
frequency. 

2.  There  is  no  slipping  between  the  elastic  and  viscoelastic  layers  at  their  interfaces. 

3.  The  major  part  of  the  damping  is  due  to  the  shearing  of  the  viscoelastic  material. 

4.  The  elastic  layers  displace  laterally  the  same  amount. 

5.  The  beam  is  simply  supported  and  vibrating  at  a  natural  frequency,  or  the  beam  is 
infinitely  long  so  that  the  end  effects  may  be  neglected. 

6.  The  normal  forces  acting  in  the  core  may  be  neglected  in  comparison  to  the  normal 
forces  in  the  elastic  material. 

7.  Only  transverse  inertia  is  considered,  the  rotatory  and  longitudinal  inertia  of  the  beam 
are  ignored. 

8.  The  face  plates  are  very  stiff  in  shear.  Therefore  shear  deformation  of  the  face  plates 
is  neglected. 
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These  assumptions  allow  for  the  closed  form  solution  to  the  viscoelastic  sandwich 
beam  problem,  but  this  has  practical  limitations  since  DiTaranto's  equation  is  in  terms  of 
only  the  longitudinal  displacements  of  the  beam  and  is  limited  to  simply  supported  end 
conditions.  Although  his  equation  is  not  directly  applicable  to  boundary  conditions  other 
than  simply  supported,  he  was  able  to  show  that  the  relationship  between  the  damping 
factor  and  the  natural  frequency  is  independent  of  the  boundary  conditions,  but  is 
dependent  upon  the  cross-sectional  geometry  and  the  physical  properties  of  the  elastic 
and  viscoelastic  materials.  Therefore  the  relationship  which  can  be  derived  for  the 
simply  supported  end  condition  will  hold  true  for  any  other  nondissipative  boundary 
conditions. 

Another  limitation  of  the  DiTaranto  equation  is  the  assumption  that  top  and  bottom 
faces  always  move  transversely  in-phase  with  one  another  ignores  possible  vibratory 
modes. 

Although  DiTaranto  did  not  complete  the  substitution,  it  can  be  shown  that  the 
differential  equation  he  solved  is  of  the  form: 


5,  +  ^3 

Rb 


K^b 


+  (fi,  +5,)f 

0 


d\ 

8x^ 


=  m 


S  d^Ui 
6~"^ 


1 

Rb  dx^dt^ 


[8] 


where  w,  =  longitudinal  displacement  of  the  top  face  layer’s  centroid 

=  /=1,3 

K,  =  E,Ai  i=l,3 
h. 

=  d 


b  =  +  +  — 

2  '  2 

K, 

R  =  ^ 
h,K, 

m  =  mass  per  unit  length 
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In  equation  [8] ,  layers  1  and  3  are  the  face  layers,  of  thickness  and  respectively, 
and  layer  2  is  the  viscoelastic  core  layer  of  thickness  h^.  and  are  the  cross-sectional 
areas  and  bending  moments  of  inertia  of  the  two  face  layers,  and  b  is  the  width  of  the 
beam. 

Use  of  the  complex  stif&iess  by  DiTaranto  in  the  equation  of  motion  for  free  vibration 
is  not  strictly  consistent  with  the  derivation  of  complex  stiffness,  which  is  based  upon 
harmonic  response.  DiTaranto  must  make  the  assumption  that  the  free  vibration  response 
is  of  the  form  =  In  the  presence  of  viscoelastic  damping  of  the  form  in 

equation  [2]  this  assumption  is  questionable. 

By  incorporating  the  complex  shear  modulus  for  the  core,  G^,  along  with  the  Young's 
moduli  for  the  two  outer  faces,  and  E2,  the  equation  contains  the  stiffiiess  of  the 
beam's  elastic  faces,  the  stiffiiess  of  the  viscoelastic  core,  and  the  damping  present  in  the 
core. 

The  resulting  derivation  considers  the  longitudinal  displacement  in  terms  of  the 
displacements  in  the  two  elastic  face  layers  along  their  centroids,  due  to  bending,  and  the 
strain  in  the  core,  due  to  shear  stress. 

Building  on  DiTaranto's  derivations  for  viscoelastic  sandwich  beam.  Mead  and 
Markus  (1969)  published  a  paper  outlining  an  equation  for  the  forced  vibration  of 
viscoelastic  sandwich  beams.  Mead  and  Markus  noted  that  DiTaranto's  sixth-order 
differential  equation  should  not  in  fact  be  applied  to  the  free  vibration  problem.  The  use 
of  complex  stiffiiess  by  DiTaranto  can  not  be  reconciled  with  the  decaying  nature  of  free 
vibration  nor  the  differential  equation  form  of  the  viscoelastic  constituitive  law.  The 
equation  however  is  applicable  to  beams  with  an  applied  harmonic  excitations,  such  as 
that  produced  by  rotating  machinery.  In  their  paper,  they  derived  an  equivalent 
formulation  of  DiTaranto's  equation  in  terms  of  transverse  displacements  instead  of  the 
longitudinal  displacements  used  by  DiTaranto.  The  derivation  was  beised  on  the  same 
assumptions  as  before.  However,  a  transverse  loading  is  applied  to  the  beam  and 
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longitudinal  and  rotatory  inertia  are  ignored.  The  solution  of  the  Mead  and  Markus 
equation  does  not  yield  the  natural  frequencies  and  damping  factors  for  viscoelastic 
sand\vich  beams  as  described  by  DiTaranto.  Instead,  the  solution  of  the  equation  yields 
special  resonant  frequencies  and  forced  modes  of  vibration.  A  closed  form  solution  for 
the  forced  vibration  problem  is  given  only  for  these  special  loading  cases  at  which  the 
modes  are  uncoupled,  in  the  sense  that  each  mode  of  vibration  is  independent  of  all  other 
modes. 

The  Mead  and  Markus  forced  vibration  equation  of  motion  is 

f  .  ,  roi 

»jl^-^j+P  =  9(x,0  [9] 

where  m  =  mass  per  unit  length 

w  =  transverse  displacement  of  the  beam  centroid 
p  =  transverse  load 

q(x,f)  =  time  dependent  transverse  external  load 

The  equation  of  motion  is  x  dependent  as  well  as  time  dependent  therefore  every  value  of 
X  represents  a  different  mode  of  vibration  for  the  beam.  The  transverse  load,/?,  contains 
the  complex  stiffness  for  the  beam  and  is  derived  by  Mead  and  Markus  as 


dS 


where  S  =  sum  of  shear  forces  in  all  three 
sandwich  layers 

The  shear  forces  in  the  two  face  plates  are  given  as 


[10] 


'  12  dx^ 

_  d^w 
^  12  dx^ 


[11] 


[12] 
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For  core  the  shear  force  is  calculated  by  multiplying  the  shear  stress  x  by  depth  d 


Sc  =  -^d 

[13] 

and 

X  =  G*y 

[14] 

where  G*  =  complex  shear  modulus  of  the  core 
y  =  shear  strain  in  the  core 

and  the  depth  (d)  is  an  approximation  which  applies  the  shear  stress  of  the 
core  over  the  thickness  of  the  core  plus  half  of  each  face  plate: 

d  =  [15] 


Mead's  Model  ACTUAL 


Figure  3  ••  Calculated  depth,  d,  compared  with 


^  x,u 


The  shear  strain  is  given  as 


_dw  ^du 
dx  dz 


[16] 
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From  the  substitution  of  equations  [10]  -  [16]  into  equation  [9]  and  the  use  of  the 
following  definitions: 


12 


Y  = 


G*(  1 


+ 


h 

(f 


\E^h^  E^h^J 


A 


\E^\+E^hj) 


[17] 

[18] 

[19] 


Along  with  the  fact  that  the  net  longitudinal  force  in  the  section  is  zero  for  transverse 
vibration,  the  sixth  order  differential  equation 


-g(i  +  J0 


m  ( 

^iDXdx^dt^ 


^  Qf2  r. 


’1  1 

■J  =  A 

^2  sq 

[20] 


is  derived  for  the  forced  vibratory  motion  of  a  beam. 

As  with  DiTaranto's  model,  the  Mead  and  Markus  equation  is  limited  in  its  practical 
application.  To  use  this  derivation,  the  special  loading  conditions  have  to  be  met. 
Although  they  derive  the  solutions  to  their  equation  for  several  different  boundary 
conditions,  the  practical  uses  are  also  limited  for  any  condition  other  than  simply 
supported  ends.  For  non-simply  supported  boundary  conditions,  an  applied  load  will 
produce  modes  of  vibration  which  may  or  may  not  be  in  phase  with  one  another  (Mead 
and  Markus,  1969).  That  means  complex  values  have  to  be  used  to  describe  these  modes 
which  greatly  increase  the  difficulty  of  the  problem. 

Lu  and  Douglas  (1974),  compared  the  predicted  harmonic  response  of  the  Mead  and 
Markus  model  with  actual  sandwich  beams.  Through  comparison  of  the  experimental 
data  with  Mead  and  Markus’s  theoretical  beam  response,  they  were  able  to  show  that 
Mead  and  Markus’s  model  predicted  the  actual  behavior  very  well  for  the  low  frequency 
response,  but  began  to  lose  accuracy  at  higher  frequencies.  More  details  of  Lu  and 
Douglas’s  paper  will  be  discussed  in  chapter  4. 
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2.2  Multiple  Variable  Models 

Although  his  paper  deals  with  multi-layered  plates  as  opposed  to  sandwich  beams,  Di 
Scuiva  (1987)  presented  a  new  displacement  model  for  layered  construction  which  offers 
some  benefit  to  this  discussion.  Recognizing  that  the  classic  plate  theory  ignores  the 
significant  effect  of  shear  deformation  on  thick  layered  plates,  Di  Scuiva  presented  a 
displacement  model  which  accounts  for  shear  deformation  in  orthotropic  plates.  The 
shear  deformation  theory  developed  prior  to  this  by  Whitney  and  Pagano  (1970)  uses  an 
arbitrary  shear  correction  factor  to  correct  for  shear  deformation  of  the  total  beam.  This 
does  not  however  take  into  account  the  effect  of  different  material  properties  in  the 
individual  layers.  In  a  multi-layered  plate,  the  deformation  pattern  is  not  linear  across  the 
entire  thickness  of  the  plate,  but  rather  it  is  only  linear  across  the  thickness  of  each  layer. 
Di  Scuiva's  model ,  known  as  the  zig-zag  theory  allows  for  each  layer  to  deform  in  a 
linear  deformation  pattern  particular  to  that  layer  and  still  allows  for  displacement 
continuity  at  the  interfaces.  Di  Scuiva's  paper  points  out  the  need  to  include  shear 
deformation  in  multilayered  material  problems  where  some  layers  have  a  much  higher 
modulus  of  elasticity  others  (Di  Scuiva,  1987). 

Averill  (1994)  adapted  Di  Scuiva's  zig-zag  model  for  plates  to  beam  analysis. 

Averill's  paper  offers  solutions  for  longitudinal,  transverse,  and  rotational  vibrations  and 
includes  the  effects  of  shear  deformation.  The  equations  offered  however  do  not  directly 
deal  with  viscoelasticity  and  therefore  do  not  consider  complex  stiffness.  For  this  reason, 
Averill's  approach  did  not  lend  itself  well  to  the  sandwich  beam  problem. 

The  paper  which  proved  to  be  the  most  beneficial  to  this  research  was  a  paper 
published  by  Bai  and  Sun  (1995).  Since  their  model  is  the  basis  for  the  finite  element 
model  developed  in  this  paper,  the  Bai  and  Sun  derivation  is  discussed  in  considerable 
detail  here.  The  effects  of  viscoelasticity  on  the  dynamic  response  of  sandwich  beams  are 
studied  by  employing  a  new  theory  which  departs  from  the  DiTaranto,  Mead  and  Markus 
models.  Bai  and  Sun  eliminate  the  assumptions  that  there  is  perfect  interface  between 
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sandwich  layers  and  that  the  top  and  bottom  faces  remain  in-plane  for  constant  transverse 
deformations.  The  core  is  modeled  as  frequency  dependent  viscoelastic  material,  while 
the  faces  are  considered  to  be  ordinary  elastic  beams  with  axial  and  bending  resistance. 
There  are  assumptions  made  about  the  conditions  at  the  interfaces  which  provide 
compatibility  and  equilibrium  requirements  to  be  satisfied  which  connect  the  behavior  of 
the  two  faces  together.  The  Bai  and  Sun  theory  also  incorporates  a  more  accurate  shear 
deformation  pattern  by  using  a  second  order  displacement  field,  T,  which  allows  the  core 
to  deform  in  a  non-linear  fashion  through  the  thickness. 


Figure  4  -  Bai  and  Sun  shear  deformation,^ 


Figure  5  -  Sandwieh  beam  Construction 
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Also  included  in  Bai  and  Sun’s  sandwich  beam  theory  was  the  consideration  of  the 
transverse,  longitudinal ,  and  rotatory  inertia  all  in  one  set  of  equations.  The  sandwich 
beam  dimensions  used  by  Bai  and  Sim  are  shown  in  Figure  5. 

Displacements  for  the  top  face  plates  are  assumed  the  same  as  for  ordinary  beams: 


u‘  =  M,  (x,t)  -  z,  — 

ox 

w'  (x,z,t)=  w^(x,t) 


[21] 

[22] 


where  u  =  longitudinal  displacement  in  the  top  face  for  any  (x,z)  location 
Ui  =  longitudinal  displacement  at  the  centroid  of  the  top  face 
Zj  =  distance  from  centroid  of  top  face  in  transverse  direction 
w=  transverse  displacement  in  the  top  face  for  any  (x,z)  location 
Wi=  transverse  displacement  at  the  centroid  of  the  top  face 


For  the  bottom  face  plates,  the  displacements  are  taken  as: 


u”  (x,z,  t)  =  (x,  t)  -  Z^  ^2  [23] 

ox 

w‘’(x,z,t)=w^(_x,t)  [24] 

where  u  =  longitudinal  displacement  in  the  bottom  face  for  any  (x,z)  location 
1/2  =  longitudinal  displacement  at  the  centroid  of  the  bottom  face 
Z2  =  distance  from  centroid  ofbottom  face  in  transverse  direction 
w’  =  transverse  displacement  in  the  bottom  face  for  any  (x,z)  location 
W2  =  transverse  displacement  at  the  centroid  of  the  bottom  face 


Displacements  for  the  core  are  assumed  as  follows 


w"  (x,  Z,  t)  =  Wq  (x,  r)  +  zO(x,  t) 


dx  ) 
z^  d^(x,t) 


2  dx  6  dx^ 


[25] 


2  dx 


[26] 
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where  u  =  longitudinal  displacement  of  core 

Uq  =  longitudinal  displacement  at  centroid  of  core 
=  transverse  displacement  of  core 
Wq  =  transverse  displacement  at  centroid  of  core 
4^=  shear  deformation  in  core 
O  =  transverse  normal  deformation  in  core 

c  =  2(1+  Vc) 

Vq  =  Poisson's  ratio  for  the  core 

In  Bai  and  Sun’s  approach  to  the  viscoelastic  sandwich  beam  problem,  an  imperfect 
adhesive  layer  is  assumed  between  the  core  and  each  of  the  faces.  To  describe  the 
interaction  at  these  interfaces,  the  adhesive  layers  are  considered  to  be  viscoelastic  layers. 
A  complex  stiffness,  1<*,  is  calculated  for  these  layers  and  the  strain  energy  contribution 
of  the  adhesive  layers  in  included  in  the  calculations  for  the  total  potential  energy  of  the 
sandwich  beam.  This  frequency  dependent  complex  stiffiiess  is  equal  to 


A*  =  k{<Si)  (l  +  y5(co)) 


[27] 


When  the  real  portion  of  the  adhesive  layer  stiffness,  ^co),  is  small,  the  adhesive  layer 
would  be  soft  and  displacements  within  this  layer  (due  to  shear  deformation)  would 
therefore  be  possible.  As  A(©)  increases  and  approaches  infinity,  there  would  be  perfect 
bonding  between  the  core  and  the  faces  just  as  all  previous  sandwich  beam  theory  has 
assumed. 

Bai  and  Sun's  research  also  suggests  that  since  the  action  within  these  adhesive  layers 
is  viscoelastic,  the  adhesive  layers  themselves  contribute  to  the  overall  damping  of  the 
sandwich  beam  system.  One  of  their  conclusions  is  that  there  may  be  an  optimal 
adhesive  layer  stiffness  that  could  be  considered  in  achieving  the  maximum  structural 
damping  for  the  system.  Because  of  its  comparatively  small  thickness,  the  effect  of 
adhesive  layer  viscoelasticity  would  not  show  up  unless  the  adhesive  layer  were 
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significantly  softer  than  the  core  itself.  An  adhesive  layer  that  is  stiffer  than  or  only 
slightly  softer  than  the  core  would  not  have  any  effect  on  the  Bai  and  Sun  calculations 
and  perfect  interface  could  be  assumed.  For  this  reason  perfect  interface  would  generally 
be  a  suitable  assumption  for  most  sandwich  beams,  unless  adhesive  layer  softness  were 
intentionally  designed  into  the  system. 

Bai  and  Sun  do  reference  a  particular  design  of  racing  boats  using  sandwiched 
materials  where  soft  adhesive  layers  increase  the  energy  absorption  of  the  system. 
However  for  the  use  of  viscoelastic  sandwich  beams  in  TMDs,  the  effort  involved  in 
designing  and  constructing  an  adhesive  layer  which  was  soft  enough  to  provide 
significant  damping  and  sturdy  enough  to  last  over  time  would  make  it  impractical.  In 
the  research  present  in  this  paper,  perfect  interface  was  assumed  in  all  calculations. 

To  "connect"  the  behavior  of  the  two  faces  and  the  core  together,  there  are  compatibility 
relationships  assumed  at  each  of  the  core-face  interfaces. 

At  the  upper  interface 

T  =  -  if)  [28] 

w^  =  \f  [29] 

where  x=  shear  stress  at  upper  interface 
At  the  lower  interface 

T  =  -  tP)  [30] 

w^  =  w^  [31] 

where  t=  shear  stress  at  lower  interface 
In  the  core 

T  =  [32] 

where  x  =  shear  stress  in  core 

G*  =  complex  shear  modulus  for  the  core 
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For  equilibrium  to  be  satisfied,  the  shear  stress  at  the  top  interface,  the  shear  stress  in  the 
core  and  the  shear  stress  at  the  bottom  interface  must  all  three  be  equal.  From  this 
equilibrium  relationship,  Bai  and  Sun  derived  the  equation: 


U^-U2  + 


/j,  +  Sw,  h2+h  dM>2 

2  dx  2  dx 


-(k*eh  +  2E,)V  +  k*—^  =  Q 

12  dx^ 


[33] 


Using  these  compatibility  equations,  Bai  and  Sim  were  able  to  rewrite  their 
displacement  field  of  eight  variables,  into  a  displacement  field  of  five  variables.  The 
details  of  this  will  be  discussed  later. 


Original  displacement  field  New  displacement  field 


Figure  7  -  Reducing  Displacement  Field 

The  governing  equations  of  motion  and  the  boundary  conditions  are  derived  via 
Hamilton's  principle,  which  requires 


h 

5  J[r-  t/J/r  =  0 

'i 


[34] 


where  T  is  the  kinetic  energy,  U  is  the  “potential  energy”,  and  and  t2  define  the  time 


interval. 
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In  the  case  of  Bai  and  Sun’s  sandwich  beam  theory,  the  kinetic  energy  includes  the 
longitudinal  inertia  of  the  top  face,  the  transverse  inertia  of  the  top  face,  the  rotatory 
inertia  of  the  top  face,  the  longitudinal  inertia  of  the  bottom  face,  the  transverse  inertia  of 
the  bottom  face,  the  rotatory  inertia  of  the  bottom  face,  the  longitudinal  inertia  of  the 
core,  and  the  transverse  inertia  of  the  core.  The  equation  for  the  kinetic  energy  of  the 
viscoelastic  sandwich  beam  then  becomes 


L 

r 

- 

(  5w,  V  2 

9  f 

2 

\pA 

("■)  +-;t 

+  P2^2 

j 

0 

12 

^5x7 

12  V  5x:  7 

+  +{wy]dzdx 

-HU 


[35] 


In  the  case  of  Bai  and  Sun’s  sandwich  beam  theory,  the  potential  energy  includes  the 
axial-extension  stiffness  of  the  top  face,  the  bending  stiffness  of  the  top  face,  the  axial- 
extension  stiffness  of  the  bottom  face,  the  bending  stifl&iess  of  the  bottom  face,  the  shear 
stiffness  of  the  core,  the  transverse  deformation  stiffiiess  of  the  core,  and  the  shear 
deformation  resistance  of  the  core.  It  is  important  to  note  that  Bai  and  Sun’s  “potential 
energy”  calculation  includes  energy  dissipation.  The  stiffness  and  the  damping  of  the 
system  are  included  together  in  one  set  of  energy  calculations  through  the  use  of  complex 
moduli.  The  equation  for  the  potential  energy  of  the  viscoelastic  sandwich  beam  then 
becomes 


U  =  U^+U^+U^^■U^+y 


where 


u,=\\ea 


5m, 


12 


dx^ 


[36] 


\dx 


[37] 
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[38] 


[39] 


[40] 


To  include  the  compatibility  requirement  of  equation  [33]  for  the  adhesive  layer,  Bai  and 
Sim  used  a  Lagrange  multiplier,  X,  which  represents  the  face-core  interface  shear  force, 
and  (j),  which  represents  the  strain  energy  within  the  adhesive  layer.  Thus  the  potential 
energy  is  augmented  by  the  expression 


L 


k 


0 


h^+h^w^  h^+hdw^ 
2  dx  2  dx 


-(k*eh  +  2E^y¥  +  k* 


12  dx^ 


[41] 


By  eliminating  Lagrange  multiplier,  X,  Bai  and  Sun  were  able  to  translate  these 
expressions  for  the  potential  and  kinetic  energy,  by  taking  variation  in  the  Hamilton's 
principle  equation  with  respect  to  ,  U2,  w  j ,  W2,  and  'F.  The  result  is  a  matrix  of  five 
governing  equations  of  motion  which  are  represented  as 


(aM)u  +  (aK)u  =f  [42] 

where  u  =  (  mj,  W2,  M2’ ^2’ 

5  M  =  the  mass  matrix  differential  operator 
5  K  =  the  stiffness  matrix  differential  operator 

f  =  (  0,  0,  0,  0  *this  is  assuming  an  applied  load,  q,  at  the  beam  center* 


It  was  also  noted  by  Bai  and  Sun  that  5  K ,  u,  and  ii  are  complex.  Some  of  the  possible 
boundary  conditions  are  derived  by  Bai  and  Sun  for  the  top  face,  bottom  face,  and  core  of 
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a  sandwich  beam  using  these  same  Hamilton's  principle  derivations.  However,  complete 
solutions  are  not  offered  for  all  boimdary  conditions  which  could  occur.  There  are 
therefore  sandwich  beam  configurations  for  which  analytical  solutions  have  not  been 
determined.  Although  Bai  and  Sun's  sandwich  beam  theory  is  the  basis  for  the  analysis 
offered  in  this  paper,  it  was  necessary  to  use  approximation  techniques  to  develop  a  more 
inclusive  sandwich  beam  analysis  tool. 
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3.  FINITE  ELEMENT  FORMULATION  FOR  DAMPED  SANDWICH  BEAMS 
3.1  Prior  Research 

The  optimal  design  of  viscoelastic  sandwich  beam  TMD's,  requires  that  an  analysis 
technique  be  available  which  could  be  used  to  find  resonant  frequencies  and  complex 
impedance  functions  for  sandwich  beams  of  any  geometric  and  material  configuration. 
Without  an  analytical  solution  offered  by  any  previous  research  which  could  satisfy  this 
requirement,  the  need  for  finite  element  modeling  became  evident. 

A  few  researchers  have  developed  finite  element  techniques  to  predict  sandwich  beam 
behavior.  Johnson,  Kienholz  and  Rogers  (1981)  developed  a  three-dimensional  model 
using  the  MSC/NASTRAN  computer  program.  They  describe  the  methods  available  for 
the  finite  analysis  of  structures  with  viscoelastic  damping  as  falling  into  three  catagories. 
The  first  is  the  complex  eigenvalue  method.  In  the  formulation  of  this  method,  the  free 
vibration  response  is  determined  using  a  complex  stiffriess  matrix,  which  is  not  frequency 
dependent.  The  mass  and  stiffness  matrices  are  taken  as  constant.  A  separate  damping 
term  can  be  added  to  account  for  viscous  damping  which  varies  with  frequency.  The 
resulting  equation  of  motion  is  solved  as  an  eigenvalue  problem  with  complex 
eigenvalues  and  eigenvectors.  There  are  two  drawbacks  to  this  approach.  First,  it  is 
difficult  to  obtain  realistic  matrices.  Very  few  viscoelastic  materials  behave  in  such  an 
accommodating  manner  (Johnson,  Kienholz,  and  Rogers,  1981).  The  second  problem  is 
computational  cost.  Complex  eigenvalue  analysis  involves  considerable  programming 
and  computing  time. 

The  next  method  type  described  by  Johnson,  Kienholz,  and  Rogers  is  the  modal  strain 
energy  method.  This  is  the  method  which  formed  the  basis  of  their  work.  By  assuming 
that  the  damped  sandvdch  beam  can  be  represented  in  terms  of  the  undamped  system  they 
were  able  to  greatly  simplify  the  viscoelastic  sandwich  beam  problem.  The  method 
involves  the  calculation  of  the  undamped  mode  shapes  of  the  sandwich  beam  with  the 
viscoelastic  material  treated  as  if  it  were  purely  elastic  with  a  real  stiffness  modulus. 
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They  then  use  the  relationship  between  the  viscoelastic  damping  factor  and  modal  strain 
energy  to  apply  modal  damping  to  the  calculations.  An  empirical  method  was  also 
included  which  allows  an  approximate  correction  for  the  frequency  dependent  material 
properties.  The  modal  strain  energy  method  was  also  incorporated  by  Soni  (1981)  in  his 
finite  element  analysis  approach  to  viscoelastic  sandwich  beams.  Soni  used  the 
MAGNA-D  computer  program  for  his  finite  element  analysis.  Both  Johnson,  Kienholz, 
and  Rogers  (1981)  and  Soni  (1981)  used  3  dimensional  brick  elements  for  their  models. 
Johnson,  Kienholz,  and  Rogers  used  an  8  node  element,  HEXA8,  for  the  core  layer,  and  a 
four  node  element,  QUAD4,  for  the  face  plate  layers.  Soni  used  isoparametric  elements 
for  his  model  which  had  8  nodes  for  the  core  layer,  and  8  to  20  nodes  for  the  face  plates 
depending  upon  the  face  plate  thickness.  In  both  papers  experimental  data  is  presented 
demonstrating  these  two  methods  on  simple  problems  with  good  results.  However  Mace 
(1992)  criticizes  these  approaches  as  being  too  complex  and  costly  to  use.  He  points  out 
the  difficulties  in  creating  the  meshes  used  and  the  complexity  of  the  analysis.  Soni  for 
instance  used  440  active  degrees  of  freedom  for  the  analysis  of  a  simple  cantilever  beam. 

The  third  finite  element  method  described  by  Johnson,  Kienholz,  and  Rogers  is  the 
direct  frequency  response  method.  This  is  the  method  used  in  this  thesis  and  was  also 
used  by  Mace  (1992).  This  method  analyzes  the  harmonically  forced  vibration  problem 
using  complex  stiffness  which  is  frequency  dependent.  The  system  damping  is  accounted 
for  in  the  complex  nature  of  the  stiffness.  Mace  directed  his  finite  element  analysis 
toward  sandwich  beams  with  very  thin  viscoelastic  layers  only.  With  the  thin  layer 
problem,  he  was  able  to  assume  the  viscoelastic  core  as  a  film  layer.  The  inertial  effects 
of  the  viscoelastic  material  are  considered  negligible  and  are  therefore  not  included  in  the 
mass  matrix.  The  stiffness  is  given  by  two  separate  matrices:  a  stiffness  matrix  for  the 
two  face  plates,  and  a  frequency  dependent,  complex  stiffhess  matrix  for  the  core.  At 
each  node,  Mace’s  model  uses  5  degrees  of  freedom,  all  of  which  represent  displacements 
in  the  two  face  plates.  Comparing  his  results  with  experimental  data  and  prior  theoretical 
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studies,  Mace  is  able  to  show  results  which  are  good,  but  only  for  studies  involving  very 
thin  core  layers.  The  finite  element  model  described  in  this  thesis  is  derived  in  much  the 
same  manner  as  the  Mace  model. 

3.2  Assumptions 

The  finite  element  model  derived  here  is  based  upon  the  assumptions  made  by  Bai  and 
Sun  (1995).  However,  some  of  Bai  and  Sun’s  assumptions  about  adhesive  layer 
behavior  were  discarded  as  unnecessary  for  the  present  study. 

Bai  and  Sun  assumed  that: 

1.  Elastic  face  layers  are  considered  to  be  ordinary  beams  with  axial  and  bending 
resistance. 

2.  The  viscoelastic  core  carries  negligible  longitudinal  stress,  but  takes  non-linear 
displacement  fields  in  both  x  and  z  directions. 

3 .  Longitudinal  displacement  discontinuity  between  the  two  face  plates  is  proportional 
to  the  shear  stress  at  the  face-core  interface,  which  is  viscoelastic. 

4.  Shear  stress  at  top  and  bottom  interfaces  equals  the  shear  stress  in  the  core  layer 

5.  Transverse  displacement  of  faces  equals  transverse  displacement  of  core  at  interfaces 

For  the  purpose  of  the  derivations  set  forth  in  this  paper,  perfect  adhesive  layer 
interaction  was  assumed  for  all  cases.  The  result  is  that  everywhere  in  Bai  and  Sun's 
calculations  where  the  adhesive  layer  stiffhess,  appears,  a  value  of  infinity  can  be 
assumed.  At  the  interfaces,  this  infinitely  stiff  adhesive  layer  transmits  the  shear  stress 
unchanged  from  the  core  to  the  face  layers. 

The  original  displacement  field  includes  eight  displacements.  Using  the  compatibility 
and  equilibrium  equations,  Bai  and  Sun  were  able  to  reduce  this  to  five  displacements. 
The  following  is  a  summary  of  how  this  was  accomplished. 


(a)  At  the  upper  interface,  the  transverse  displacement,  w,  of  the  core  is  equal  to  the 
transverse  displacement  of  the  top  face.  Substituting  equations  [22]  and  [26]  into 
equation  [29]  with 
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yields  the  equation 


w,  =  Wfl  +—<!>- 
1  0  2 


8  dx 


[43] 


(b)  At  the  lower  interface,  the  transverse  displacement,  w,  of  the  core  is  equal  the 
transverse  displacement  of  the  bottom  face.  Substituting  equations  [24]  and  [26]  into 
equation  [31]  with 


yields  the  equation 


8  dx 


(c)  Solving  equations  [43]  and  [44]  for  and  O  yields 


and 


>^0 


w,  +  w, 

— - -  + - 

2  S  dx 


<D  = 


w,  -  W2 


h 


[44] 


[45] 


[46] 


In  the  remainder  of  the  discussion,  a  simplified  partial  differential  equation  notation  will 
be  used,  to  reduce  the  length  of  the  equations 
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d{.) 

etc 

becomes 

( •  )»x 

a^(.) 

dx^ 

becomes 

( •  )’XX 

5^(.) 

dx^ 

becomes 

( •  )»xxx 

Thus,  for  instance  Equation  [45]  becomes 


w,  +  Wj 
2 


(d)  At  the  upper  interface,  the  longitudinal  displacement,  u,  of  the  core  equals  the 
longitudinal  displacement  of  the  top  face.  Substituting  equations  [21]  and  [25]  into 
equation  [28]  with 


yields  the  equation 


k  ^  ^  ^  .Tr 


[47] 


(e)  At  the  lower  interface,  the  longitudinal  displacement,  u,  of  the  core  equals  the 
longitudinal  displacement  of  the  bottom  face.  Substituting  equations  [23]  and  [25]  into 
equation  [30]  with 


u 


c 


yields  the  equation 


29 


^2  2  ^2  5JC  ^0  2  ^0  5a:  )  g  ®?;c  ^g  ^5 


[48] 


(f)  Solving  equations  [47]  and  [48]  allows  Woto  be  expressed  in  terms  of  the  face  plate 
displacements  only  as 
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The  displacement  field  of  eight  variables  can  now  be  represented  by  five  variables.  To 
reduce  the  size  of  the  problem  even  further  the  last  compatibility  expression,  equation 
[33]  is  used.  When  the  adhesive  layer  stiffiiess,  k*,  is  taken  to  be  infinite,  this  expression 
can  be  re-written  as 
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[50] 


Using  this  expression,  it  is  possible  to  eliminate  from  the  displacement  field  the 
variable  T.  Initial  attempts  to  include  T  as  a  variable  in  the  finite  element  formulation 
were  unsuccessful,  resulting  in  badly  behaved  finite  element  models.  However,  insight 
was  gained  fi'om  those  initial  attempts  that  permitted  an  approximate  model  to  be 
constructed  successfully. 

The  first  attempted  element  derivation  simply  replaced  in  the  variational 
formulation  with  the  expression 


vp  _  l^vp  . 
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[51] 
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Then 'T  was  approximated  using  linear  shape  functions  as 

'J'(X,0  =  [l-2]'P,+^'P;  [52] 

within  the  element.  This  approximation  led  to  a  badly  behaved  element  whose  stiffness 
appeared  to  increase  as  the  number  of  elements  were  increased. 

A  second  attempt  to  include  explicitly  as  a  nodal  variable  was  likewise 
imsuccessfiil,  but  suggested  the  final  form  of  the  model.  Equation  [50]  is  a  differential 
equation  for  'F(x,t)  in  terms  of  the  face  plate  displacements  and  their  derivatives. 
Assuming  harmonic  response  of  u^,  u^,  Wj,  W2,  and  T  allows  equation  [50]  to  be 
rewritten  as  an  ordinary  differential  equation  describing  the  spatial  variation  of 'E  at  any 
given  time,  t.  The  solution  of  this  differential  equation  is  given  as  the  sum  of  the 
homogeneous  solution  'E;,  and  the  particular  solution 
Taken  as  an  ordinary  differential  equation,  equation  [50]  is  quite  easy  to  handle,  as  the 
coefficients  are  constant.  The  homogeneous  solution  is  of  the  form 

iHfjr  - 

'E^^Cx)  =  +026^*'  [53] 


which  can  also  be  written  as 


'E*  (x)  =  Aj  cosh 


X  +  .^2  cosh 


[54] 


The  homogeneous  solution  can  be  used  to  construct  a  particular  solution,  by  the  method 
of  variation  of  parameters,  for  given  U2,  Wj,  and  W2  functions.  However  a  more 
straight  forward  approach  is  to  recognize  that  the  homogeneous  solution  contains  all 


information  needed  to  construct  the  nodal  equations,  since  it  contains  the  boundai^' 
information.  Thus,  it  was  decided  to  construct  a  new  set  of  shape  functions  for 
using  the  homogeneous  solution  of  equation  [54].  This  leads  to  the  new  equations 
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This  shape  flmction  for  also  proved  to  be  unsatisfactory,  and  in  fact  led  to  singular 

equations.  Investigation  into  the  variation  of  the  shape  functions  of  equation  [55] 

revealed  that  the  influence  of  boundary  conditions  upon  the  function  extends  only 

a  very  short  distance  into  the  beam  for  most  realistic  values  of  — .  Figure  8  illustrates 

L 

this  point.  On  this  figure,  the  value  of  the  left  hand  shape  function  as  a  function  of  —  is 

L 

h 

plotted  for  various  values  of  — .  It  is  seen  firom  this  figure  that  the  shape  function 

L 

effectively  reaches  zero  for  y  <  0.1  and  —  >  about  0.1.  Thus  for  the  relatively  thin 

L  L 

sandwiches  likely  to  be  used  in  the  development  of  a  viscoelastically  damped  sandwich 
beams,  boundary  conditions  do  not  propagate  very  far  into  the  beam.  Similarly,  any  local 
variation  of  T  within  the  beam  should  not  propagate.  This  suggests  that  the  second 
derivative  term  in  equation  [50]  is  of  minor  importance,  and  we  may  approximately  state 
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[56] 


This  approximation  should  be  satisfactory  provided  the  ratio 


A 

L 


is  sufficiently  small. 


As  the  beam  is  subdivided  into  small  beam  elements  the  ratio  of  core  thickness  to  beam 
h 

element  length,  — ,  increases  which  suggests  the  number  of  elements  that  can  effectively 
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be  used  for  a  given  beam,  before  this  assumption  becomes  invalid,  has  a  limit.  For 
example,  if  a  beam  of  length  50  inches  has  a  core  thickness  of  .5  inches,  the  ratio  is 

.01  which  is  small.  If  during  the  finite  element  process  the  beam  were  discretized  into  50 
h 

elements,  the  ratio  —  would  go  up  to  .5  for  each  of  these  beam  elements.  This  suggests 
L 

that  the  homogeneous  portion  of  the  shear  deformation  would  become  important  for  the 
entire  length  of  the  beam  element.  Thus,  two  important  questions  to  be  addressed  in  the 
analysis  is  whether  sufficiently  accurate  frequency  response  information  can  be  obtained 
without  resorting  to  this  level  of  discretization,  and  whether  the  approximation  does  in 
fact  cause  the  element  to  perform  badly  for  fine  discretization  meshes. 

With  the  approximation  introduced  by  equation  [56]  it  is  possible  to  eliminate  T  from  the 
model.  The  displacement  field  is  now  expressed  in  terms  of  four  displacements:  the 
transverse  and  longitudinal  displacements  of  the  top  face  plate  (Wj  and  Wi)  and  the 
transverse  and  longitudinal  displacements  of  the  bottom  face  plate  {u2  and  W2).  Since  the 
compatibility  equation  [33]  has  been  explicitly  accounted  for,  it  is  further  possible  to 
eliminate  the  term  ^  from  the  total  potential  energy. 

3.3  Assuming  Shape  Functions 

In  their  equations  for  potential  and  kinetic  energy,  Bai  and  Sun  define  the  energy 
equations  in  terms  of  the  displacements  and  the  partial  derivatives  of  the  displacements 
with  respect  to  x.  Because  these  displacements  are  both  functions  of  x  and  t,  the  solution 
of  these  equations  requires  that  these  displacements  be  separated  into  two  distinct 
displacements-  one  a  function  of  x  and  one  a  function  of  t. 

In  general  the  separation  of  variables  takes  the  form 


M,(x,0  =  a,(x)M,(0 


[57] 


>^1(^.0  =  Pi  (^M(0 


[58] 
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[59] 

W2(x,0  =  P2WW2(0 

[60] 

In  the  case  of  the  model  developed  for  harmonic  response,  these  equations  can  be  reduced 
to  the  form 


u^{x,t)  -  a,(x)e“-" 

[61] 

Wi(x,0  =  PiCxK-'' 

[62] 

M2(x,0  =  ct2(x)e“-'' 

[63] 

W2(x,0  =  p2We“-'' 

[64] 

It  is  thus  possible  to  solve  the  energy  equations  as  functions  of  x,  for  a  given  co.  Once  the 
functions  a ,  (x) ,  p,  (x) ,  ttj  (x) ,  and  P2  (x)  are  specified  in  terms  of  shape  functions  and 
nodal  values  this  solution  yields  the  mass  and  stiffhess  matrices  for  the  sandwich  beam 
system.  The  shape  functions  are  chosen  in  the  usual  manner. 

In  the  Bai  and  Sim  potential  and  kinetic  energy  equations,  the  displacements  are 
defined  in  terms  of  partial  derivatives  of  the  displacements.  The  displacement  field  has 
been  reduced  to  four  variables.  The  two  longitudinal  displacements,  u,  show  up  in  the 
energy  equation  as  first  order  derivatives.  The  shape  functions  chosen  to  represent  them 
are  the  standard  Q  shape  functions  used  for  truss  elements  and  are  shown  in  Figure  9. 

The  two  transverse  displacements,  w,  used  in  this  derivation  show  up  in  the  energy 
equation  as  second  order  derivatives.  The  shape  fimctions  chosen  to  represent  them  are 
the  standard  beam  shape  fimctions,  shown  in  Figure  10.  Third  order  derivatives  for  the 
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transverse  displacements  appear  in  the  kinetic  energy  equations  once  all  substitutions 
have  been  made.  This  however  shows  up  only  as  part  of  the  term  for  the  longitudinal 
inertia  of  the  core.  Since  the  core  density  is  usually  significantly  smaller  than  the  face 
density,  this  inertia  term  is  small.  The  third  derivative  terms  therefore  have  little  effect 
on  the  behavior  of  the  beam  system  and  were  not  considered  when  determining  shape 
fimctions. 

When  the  interpolation  fimctions  are  placed  into  the  energy  equations,  it  is  then 
possible  to  carry-out  the  integration  of  the  energy  equations  and  form  the  mass  and 
stiffness  matrices.  With  the  inclusion  of  the  interpolation  fimctions,  the  displacement 
field  becomes  a  vector  of  six  displacements  at  each  beam  node. 


Figure  11  -  Displacement  Vector 

With  twelve  displacements  to  be  calculated,  the  element  mass  and  stiffiiess  matrices 
become  twelve  by  twelve  matrices.  It  is  impractical  to  include  all  the  details  of  the  mass 
and  stiffiiess  matrix  formulations  in  this  thesis  since  the  work  involved  in  assembling 
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these  twelve  by  twelve  matrices  is  long  and  tedious.  It  should  be  noted  that  since  the  top 
and  bottom  face  plates  are  taken  to  behave  as  ordinary  elastic  beams,  portions  of  these 
matrices  correspond  to  standard  beam  matrices  where  the  displacements  {u,  w,  and  w ' ) 
are  considered.  The  complete  matrices  are  provided  in  Appendix  1 .  The  assembled 
matrices  are  also  included  as  part  of  the  FORTRAN  program  located  in  Appendix  2. 

Because  of  the  uniqueness  of  the  viscoelastic  sandwich  beam  system,  pre-developed 
finite  element  codes  typically  do  not  address  the  structural  analysis  of  such  systems.  To 
perform  the  discretization  and  assembly  process  of  the  finite  element  method,  it  was 
therefore  necessary  to  write  a  computer  program  especially  for  the  sandwich  beam 
problem.  This  was  accomplished  using  the  FORTRAN  programming  language. 

3.4  FORTRAN  Program 

The  FORTRAN  program  was  written  with  the  mass  and  stiffiiess  matrices  included  in 
the  programming.  These  matrices  are  provided  with  unknown  variables  allowing  the 
program  to  model  any  sandwich  beam  configuration.  No  attempt  was  made  at  this  stage 
of  element  development  to  develop  a  general  code,  but  the  element  developed  should  be 
easily  incorporated  into  such  a  program. 

Once  the  sandwich  beam  information  is  provided,  the  mass  and  stiffness  matrices  for 
each  element  are  calculated  and  stored.  The  program  then  must  assemble  each  of  the 
individual  element  matrices  into  one  a  global  mass  and  stiffness  matrix  for  the  entire 
beam.  By  defining  the  displacement  field  so  that  the  displacements  for  each  end  node 
are  kept  together,  it  was  possible  to  add  these  matrices  in  a  convenient  manner.  The 
displacements  for  node  j  for  one  element  are  the  same  as  the  displacements  for  node  /  for 
the  next  element  and  so  forth. 

3.5  Equation  of  Motion 

Once  mass  and  stiffiiess  matrices  were  derived,  the  equation  of  motion  could  be 
determined.  Using  the  direct  frequency  response  method  the  form  of  the  equation  of 
motion  for  a  beam  element  will  be 
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Mu  +  K*u  =  q(/)  [65] 

where  M  =  mass  matrix 

u  =  nodal  displacement  vector 
ii  =  nodal  acceleration  vector 
K*  =  complex  stiffiiess  matrix,  dependent  upon  Q 
q(r)  =  nodal  force  vector 

If  a  harmonic  load  is  applied 

q(0  =  [66] 

where  Q  is  the  forcing  frequency.  Both  amplitude  and  mode  of  vibration  are  dependent 
on  the  frequency.  It  can  be  assumed  that  the  response  due  to  this  applied  harmonic  load 
will  also  be  harmonic  and  at  the  same  frequency,  Q.  Then 

u  = 

ii  = 

The  equation  of  motion  then  becomes 

(K*  -  Q^M)  U  =  Fq 

where  U  =  (  u^,  wjy,  w’^,  M2,-,  1^2/’  ^2/’  "1/’  ^1/’  ^’ip  "2/’  ^2p  ^2/ 

M  =  the  mass  matrix 

K  =  the  stiffness  matrix 

Fo  =  ( 0,  Q,  0, 0, 0, 0, 0, 0, 0,  0, 0, 0,  0, 0  )T 

This  equation  of  motion  then  describes  the  response  of  the  sandwich  beam  with  an 

applied  harmonic  load  and  likewise  the  response  of  one  beam  element. 

The  FORTRAN  program,  uses  an  iterative  process  to  combine  the  mass  and  stiffness 

matrices  over  a  range  of  forcing  frequencies.  Because  the  stiffhess  matrix  is  frequency 

dependent,  the  matrix  must  be  recalculated  for  each  new  frequency. 

3.6  FORTRAN  Solution  Process 

Using  inputted  information  about  a  sandwich  beam’s  geometry  and  material 
properties,  the  FORTRAN  program  is  able  to  solve  the  sandwich  beam  problem.  The 
first  step  in  this  solution  process  is  to  calculate  the  mass  matrix.  The  mass  matrix,  M,  is 
not  frequency  dependent  and  therefore  constant.  Next  the  program  begins  an  iterative 


[67] 

[68] 

[69] 
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process  for  each  frequency  to  be  considered.  For  a  given  frequency,  Q,  the  core 
properties  are  calculated  using  the  appropriate  material-frequency  function.  Using  these 
core  properties,  the  complex  stiffness  matrix,  K*,  is  updated.  The  matrix  K*  -  Q  M  is 
calculated.  Using  the  IMSL  subroutine  LINCG  to  invert  the  matrix,  equation  [69]  is 
solved  for  the  nodal  displacements  at  that  frequency. 
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4.  NUMERICAL  STUDIES 

To  test  the  validity  of  results  from  the  FORTRAN  program,  the  system  was  used  to 
analyze  two  previously  investigated  sandwich  beam  problems  and  results  were  compared. 
Also,  by  using  Ansys,  a  computer  structural  analysis  program,  a  simplified  model  was 
constructed  of  a  sandwich  beam  in  which  the  complex  nature  of  the  core  was  ignored. 
Ansys,  which  is  also  a  finite  element  based  approach,  was  able  to  determine  the 
frequencies  of  the  first  five  modes  of  free  vibration.  Results  were  then  compared  by 
using  the  FORTRAN  program  to  find  the  undamped  vibration  response  for  the  same 
beam. 

4.1  Test  Cases 

The  test  cases  were  provided  in  a  paper  written  by  Lu  and  Douglas  (1974).  The 
purpose  of  their  paper  was  to  compare  Mead  and  Markus’s  sandwich  beam  analysis  to 
actual  test  cases.  Bai  and  Sun  (1995)  then  used  the  results  found  by  Lu  and  Douglas  to 
compare  their  model  to  both  Mead  and  Markus’s  model  and  actual  test  data.  The  Mead 
and  Markus  analysis  provided  very  accurate  results  for  comparisons  on  a  thin  core 
layered  sandwich  beam.  For  comparisons  on  a  thick  core  layered  beam  their  results  were 
again  accurate  for  the  first  and  second  modes  of  vibration.  After  that  the  results  began  to 
skew  from  the  actual  test  data  and  became  worse  as  each  higher  mode  was  reached.  Bai 
and  Sim’s  model  however  remained  consistently  accurate  even  at  the  higher  modes  of 
excitation.  The  two  test  cases  were  based  on  a  sandvvich  beam  with  the  characteristics  as 
shown  in  Table  1 .  The  core  material  used  in  specimen  1  was  an  acrylic  based  material 
which  had  a  mass  density  of  .000103  Ib-secVin'*.  The  frequency  dependent  core  complex 
shear  modulus  was  calculated  at  a  temperature  of  124°  F  to  be 

G*  =  G'+  JG"  [70] 

where  G’ =  0.579 In/ +  1.136  Ibf/in^ 

G"  =  0.6011n/  + 1.144  Ibfrin^ 

/  =  frequency  in  a  range  of  50  to  3000  Hz 
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The  core  material  used  in  specimen  2  was  neoprene  with  a  mass  density  of  .0001 1 5  Ib- 
sec^/in"*.  The  frequency  dependent  core  complex  shear  modulus  was  calculated  at  a 
temperature  of  77°  F  to  be 


G*  =  G(l  +  jb) 


[71] 


where  G  =  0.0261n/  + 4.754  IbfW 
5  -0.1 70 In/ +  2.705  IbfW 
/  =  frequency  in  a  range  of  50  to  10000  Hz 

The  analysis  was  done  over  a  range  of  frequencies  for  an  applied  load  located  at  the 
beam  center.  To  show  results,  this  range  of  frequencies  was  plotted  against  the  resulting 
driving  point  mechanical  impedance  of  the  sandwich  beam  system.  Mechanical 
impedance  is  defined  as  the  force  of  the  applied  loading  over  the  velocity  of  the  system 
in  the  direction  of  the  load.  For  the  calculation  of  impedance  for  the  test  cases  the  force 
is  given  as  a  vector  with  only  one  non-zero  corresponding  to  a  concentrated  nodal  forcing 
of  unit  magnitude. 


Specimen  1 

Specimen  2 

face  material: 

Steel 

Steel 

core  material: 

Aciylic  base 

Neoprene 

core  thickness: 

.004  inches 

.25  inches 

top  face  thickness: 

,25  inches 

.25  inches 

bottom  face  thickness: 

.25  inches 

.25  inches 

beam  width: 

1 .0  inches 

2.0  inches 

beam  length: 

24.1875  inches 

18.125  inches 

boundary  conditions: 

free-free  ends 

free-free  ends 

Table  1  -  Test  Case  Characteristics 

Using  the  given  equation  of  motion  from  equation  [69],  the  mechanical  impedance  at 
degree  of  freedom  i  under  forcing  corresponding  to  degree  of  freedom  k  can  be  defined  as 
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ImPi,  =  — ^  [72] 

jDU,. 

The  amplitude  of  this  impedance  is 

|lmp^|=  [73] 

where  U,  =  complex  conjugate  of  U, 

Using  the  FORTRAN  program  to  analyze  the  sandwich  beam  described  above  the 
mechanical  impedance  was  calculated  at  several  different  levels  of  beam  discretization. 
Plotted  in  Figures  12  and  14  are  the  results.  Figure  13  shows  Lu  and  Douglas’s 
comparison  of  experimental  results  with  Mead  and  Markus’s  analysis. 

The  specimen  1  sandwich  beam  had  a  very  thin  core  of  viscoelastic  material.  Since 
Bai  and  Sim  were  primarily  interested  in  thick  layer  sandwich  beams,  they  did  not  use 
this  data  to  compare  their  model.  Analysis  of  the  beam  by  the  FORTRAN  program  was 
done  using  an  8  element  model  (taking  advantage  of  beam  symmetry  the  beam  was  only 
modeled  to  the  midpoint  for  all  meshes  discussed  here).  Results  are  plotted  in  Figure  12 
where  the  red  line  represents  the  theoretical  results  and  the  boxes  show  the 
experimentally  obtained  data.  The  results  show  the  FORTRAN  model  has  a  great  deal  of 
accuracy  in  matching  the  experimentally  obtained  data,  especially  at  the  lower 
frequencies. 

The  specimen  2  sandwich  beam  had  a  thick  core  of  a  very  soft  viscoelastic  material. 
The  results  of  Lu  and  Douglas’s  experiments  with  specimen  2  are  plotted  in  Figure  13. 
This  plot,  showing  the  comparison  of  the  experimentally  obtained  data  and  Mead  and 
Markus’s  theory  based  results,  demonstrates  how  Mead  and  Markus’s  model  begins  to 
skew  at  higher  frequencies.  Analysis  of  the  beam  was  performed  by  the  FORTRAN 
model  and  the  results  plotted  in  Figure  14  for  a  2  element  mesh  (represented  as  a  blue 
line)  and  an  8  element  mesh  (represented  as  a  red  line).  The  analysis  provided  accurate 
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Figure  12  -  Specimen  1  Test  Data 


Line  -  8  element  mesh 
Boxes  -  Experimental  Data 
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Test  Data  &  Program  Results 


Figure  14  -  Specimen  2  Test  Data 


Blue  Line  -  2  element  mesh 
Red  Line  -  8  element  mesh 
Boxes  -  Experimental  Data 
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results  for  lower  frequency  response,  even  with  the  very  crude  2  element  model.  As  the 
number  of  elements  were  increased,  the  model  began  to  demonstrate  accuracy  even  at  the 
higher  frequencies.  In  fact,  tvithin  the  range  of  frequencies  plotted,  a  4  element  model 
provided  results  that  were  very  close  to  the  8  element  results. 

For  additional  results,  an  Ansys  computer  model  was  generated  using  the  program  in 
Appendix  3.  This  is  a  model  of  a  beam  similar  to  specimen  2  in  geometry,  however,  in 
this  analysis  the  damping  factor  was  considered  to  be  zero.  In  addition,  the  frequency 
dependence  of  the  real  part  to  the  core  moduli  was  ignored.  The  core  stiffiiess  was 
calculated  as  a  constant  equal  to  the  specimen  2  core  stiffness  at  a  frequency  of  zero.  The 
beam  mesh  used  is  shown  in  Figure  15.  Ansys,  considering  the  beam  as  a  simply 
supported  undamped  system,  calculated  the  first  five  modes  of  free  vibration.  As  a 
comparison  to  these  results,  the  FORTRAN  model  was  used  to  analyze  a  similar  beam 
problem.  The  damping  factor  was  set  to  zero  as  with  the  Ansys  model,  however  the 
frequency  dependence  of  the  real  portion  to  the  core  stiffness  remained.  Analysis  of  this 
beam  was  performed  as  a  harmonically  loaded,  simply  supported  beam.  Although  the 
differences  in  the  analysis  does  not  allow  for  an  exact  comparison  of  the  Ansys  and 
FORTRAN  program  results,  a  rough  comparison  demonstrates  that  the  FORTRAN 
program  provided  modes  of  vibration  which  are  analogous  to  the  Ansys  results.  These 
results  are  included  in  Table  2. 


ANSYS  MODEL 

FORTRAN  MODEL 

MODE 

FREQUENCY  (Hz) 

FREQUENCY(Hz) 

1 

65.75981 

66 

2 

418.15008 

408 

3 

702.48751 

667 

4 

765.98113 

759 

5 

994.21499 

983 

Table  2  -  Ansys  results 
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Figure  15  -  Ansys  model 
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5.  CONCLUSIONS 

The  theoretical  studies  provided  have  shown  that  the  finite  element  model  developed 
and  described  in  this  thesis  appears  to  be  an  effective  analysis  tool.  Both  thin  and  thick 
layer  beams  were  analyzed  with  equally  promising  results.  Although  data  was  not 
available  to  compare  with  this  model’s  higher  frequency  results  (only  data  for  frequencies 
up  to  3000  Hz  were  available),  the  results  that  are  available  show  a  definite  pattern  of 
convergence  towards  accurate  results.  The  use  of  this  model  for  lower  frequency 
analysis,  which  would  constitute  its  predominant  use,  has  been  shown  to  provide  good 
results. 

Because  of  the  assumptions  about  shear  deformation  behavior,  it  was  assumed  that  it 
is  necessary  to  limit  the  number  of  elements  used  in  any  mesh,  so  that  the  core  thickness 
to  beam  length  ratio  does  not  become  large.  The  theoretical  studies  provided  encouraging 
information  about  this  assumption  as  well.  The  model  appears  to  converge  very  quickly 
for  even  crude  meshes.  This  is  important  since  it  indicates  that  the  use  of  this  program 
may  never  require  a  large  number  of  elements  to  achieve  accurate  results.  Also  important 
is  that  as  the  number  of  elements  is  increased,  the  results  continue  to  converge  towards 
accuracy.  A  point  was  never  reached  at  which  the  number  of  elements  in  the  mesh  began 
to  diverge  away  from  the  correct  results.  This  indicates  that  even  at  higher  core  thickness 
to  beam  length  ratios,  the  assumed  shear  deformation  may  provide  good  results.  Further 
investigation  would  need  to  be  performed  to  find  if  there  were  a  limit  at  which  the 
model’s  assumptions  became  invalid  making  the  model  unsuitable.  With  results 
converging  as  quickly  as  they  do,  a  large  number  of  beam  elements  would  appear  only  to 
be  necessary  in  very  rare  instances  such  as  the  analysis  of  very  high  frequency  responses. 

Although  a  more  in-depth  study  of  this  model  will  be  necessary  before  its  benefits  are 
truly  known,  the  future  application  of  this  model  is  favorable.  The  FORTRAN  model 
provided  can  be  used  as  the  basis  for  building  future  programs  designed  to  perform  many 


types  of  analysis.  The  use  of  this  model  should  be  particularly  important  to  the 
development  of  viscoelastic  sandwich  beams  for  use  as  tuned  mass  dampers. 
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APPENDIX  1 
Mass  and  Stiffiiess  Matrices 


The  complex  stiffness  matrix ,  K*  =  K,  +  K2  +  K3  +  K4 
and 

the  mass  matrix ,  M  =  Mj  +  M2  +  M3 


where  Kj,  K2,  K3,  K4,  Mj,  M2,  and  M3  are  shown  on  the  following  pages 
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APPENDIX  2 

FORTRAN  Program 


PROGRAM  SANDWICH 

INTEGER  BIT,  ROW,  COL,  ICONNO(25) 

REAL  Mass(  12,1 2),lth,Mmatrix(606,606),Imped 
REALK,L,ST1(12,12),M1(12,12),M2(12,12),M3(12,12) 

COMPLEX  Ec,  C,  D,  E,  F,  G,  H,  Gc 
COMPLEX  Kmatrix(606,606),  Sil,Si2 

COMPLEX  STIFFNESS(12, 12),ST2(12,12),ST3(12,12),ST4(12,12) 

COMPLEX  COMBINED(606,606),  CINV(606,606) 

COMPLEX  Amatrix(606),Aconj 

PRINr(lX,A)','  V  V  V  ' 

PRINr(lX,A)',’ V  V  ’ 

PRINT'(29X,A)','SANDWICH  BEAM  ANALYSIS' 

PRINT  ’ 

PRJNT'(14X,A)','Finite  Element  Model  for  3-Layered  Sandwich  Beams' 

PRINT  '(1X,A)',' ',' ',' ' 

PRINT  *,'PLEASE  ENTER  THE  NUMBER  OF  ELEMENTS  TO  BE  USED' 
READ  *,NEL 
PRINT  *,' ' 

PRINT*, 'PLEASE  ENTER  VALUE  FOR  THE  BEAM  LENGTH  (in  INCHES):' 
READ  *,LTH 

PRINT  *,'BEAM  LENGTH  IS  ',LTH,'  inches' 

PRINT  *,' ' 

PRINT*  ,'PLEASE  ENTER  VALUE  FOR  THE  TOP  FACE  PLATE’, 

&'  THICKNESS  (in  INCHES):' 

READ  *,H1 

PRINT  *,'TOP  FACE  THICKNESS  IS  ',H1,'  inches’ 

PRINT  *,' ' 

PRINT*  ,’PLEASE  ENTER  VALUE  FOR  THE  BOTTOM  FACE  PLATE', 

&’  THICKNESS  (in  INCHES):' 

READ  *,H2 

PRINT  *,’BOTTOM  FACE  THICKNESS  IS  ',H2,'  inches' 

PRINT  *,' ' 

PRINT*,'PLEASE  ENTER  VALUE  FOR  THE  CORE  THICKNESS  (in  INCHES):’ 
READ  *,HH 

PRINT  *,'CORE  THICKNESS  IS  ',HH,'  inches’ 

PRINT  *,'  ’ 

PRINT*,'PLEASE  ENTER  VALUE  FOR  THE  BEAM  WIDTH  (in  INCHES):’ 


READ*,WW 

PRINT  *,’BEAM  WIDTH  IS  ',WW,’  inches’ 

PRINT  *; ' 

*  Beam  Length  (divided  into  elements) 

1th  =  lth/REAL(nel) 

PRINT  *, 'INPUT  NUMBER  OF  CONSTRAINTS  ON  TOTAL  BEAM’ 
READ  *,NCONSTR 
DO  25,  ROW  =  LNCONSTR 

PRINT  *, ’INPUT  #’,ROW,’  ROW  TO  BE  CONSTRAINED’ 

READ  *,  ICONNO(ROW) 

25  CONTINUE 

PRINT*,’INPUT  MAGNITUDE  OF  POINT  LOAD’ 

READ  *,Fo 

PRINT*, ’LOAD  MAGNITUDE  IS  ’,Fo 

PRINT*,’INPUT  ROW  NUMBER  WHERE  POINT  LOAD  IS  APPLIED’ 
READ  *,NMS 

PRINT  *,’TOP  FACE  PLATE  MATERIAL  IS  (Choose  one):’ 

PRINT  ’(35X,A)’,  ’1.  STEEL’ 

PRINT  ’(35X,A)’,  ’2.  OTHER’ 

READ  *,BIT 
IF  (BIT  .EQ.  1)  THEN 
CALL  STEEL(rl,el) 

ELSE 

CALL  OTHER(rl,el) 

END  IF 

PRINT  *,’BOTTOM  FACE  PLATE  MATERIAL  IS  (Choose  one):’ 
PRINT  ’(35X,A)’,  ’1.  STEEL’ 

PRINT  ’(35X,A)’,  ’2.  OTHER’ 

READ  *,BIT 
IF  (BIT  .EQ.  1)  THEN 
CALL  STEEL(r2,e2) 

ELSE 

CALL  OTHER(r2,e2) 

END  IF 

PRINT  *,’CORE  MATERIAL  IS  (Choose  one):’ 

PRINT  ’(35X,A)’,  ’1.  NEOPRENE’ 

PRINT  ’(35X,A)’,  ’2.  OTHER’ 

READ  *,NET 

IF  (NET  .EQ.  1)  THEN 

EVAL  =  2.98 

rr  =  .000115 

ELSE 

PRINT* ,’INPUT  POISSONS  RATIO:’ 

READ  *,VC 

PRINT*  ,’INPUT  DENSITY:’ 


A2-3 


READ  ♦^RR 
EVAL  =  2.0*(1.0+VC) 
END  IF 


*  mass  density  of  steel  multiplied  by  width  of  beam 

rl  =rl*WW 
r2  =  r2*WW 

*  mass  density  of  core  (x  width) 

rr  =  rr*WW 


El  =E1*WW 
E2  =  E2*WW 
A  =  (El  *  hi) /1th 
B  =  (E2  *  h2)  /  1th 
K  =  (A*  (hi /1th)**  2) 

L  =  ((A*(hl  **  2))/ 2.0 /1th) 
P  =  ((A*(hl  **2))/3.0) 

R  =  ((A  *  (hi  **  2))  /  6.0) 

S  =  (B  *  (h2  /  1th)  **  2) 

T  =  ((B*(h2**2))/2.0/lth) 
W  =  ((B*(h2**2))/3.0) 

Y  =  ((B  *  (h2  **  2))  /  6.0) 
ALPHA=0.0 


*23456789012345678902234567890323456789042345678905234567890623456789072 

*  Input  Mass  Matrix 

CALL  Massl(hh,hl,h2,lth,rr,eval,Ml) 

CALL  Mass2(hh,rr,eval,hl  ,h2,lth,Alpha,M2) 

CALL  Mass3(rl,r2,hl,h2,lth,M3) 

DO  53,ROW=l,12 
DO  54,  COL=l,12 

Mass(ROW,COL)=Ml(ROW,COL)+M2(ROW,COL)+M3(ROW,COL) 

54  continue 
53  continue 

*USE  THIS  COMMAND  TO  PRINT  THE  MASS  MATRIX  AND  ITS  COMPONENTS 

*  CALL  WRRRN  ('Mass-wc',  12, 12,  Ml,  12, 0) 

*  CALL  WRRRN  ('Mass-uc',  12, 12,  M2, 12, 0) 

*  CALL  WRRRN  (’Mass-w  &  u',  12, 12,  M3, 12,  0) 

*  CALL  WRRRN  ('Mass’,  12, 12,  Mass,  12, 0) 

DO  3000,  BIT=1,183 


IF  (BIT  .LE.  50)  THEN 


FREQ=BIT*2.0 
ELSE  IF  (BIT  .LE.  75)  THEN 
FREQ=(BIT-25.0)*4.0 
ELSE  IF  (BIT  .LE.  85)  THEN 
FREQ=(BIT-55.0)*10.0 
ELSE  IF  (BIT  .LE.  120)  THEN 
FREQ=(BIT-70.0)*20.0 
ELSE  IF  (BIT  .LE.  145)  THEN 
FREQ=(BIT-95.0)*40.0 
ELSE  IF  (BIT  .LE.  183)  THEN 
FREQ=(BIT-120.0)*80.0 
END  IF 

IF  (NET  .EQ.  1)  THEN 

CALL  NEOPRENE(FREQ,RR,EC,VC) 

ELSE 

CALL  OCORE(RR,EC) 

END  IF 

IF  (BIT  .EQ.  1)  THEN 

PRINT  DENSITY  YOUNGS  MODULUS', 

&'  POISSONS  RATIO' 

PRINT  '(30X,A)','REAL  IMAGINARY' 

PRINT  21, 'TOP  FACE  ',Rl/ww,El/ww 
PRINT  20, 'CORE  ',RR,EC,'  ',VC 

PRINT  21  ,'BOTTOM  FACE  ',R2/ww,E2/ww 

20  FORMAT  (A16,E9.3,2  F13.3,A9,F4.2) 

21  FORMAT(A16,E9.3,F13.3) 

END  IF 

Ec=Ec*ww 

C  =  (13.0  *Ec*  1th)/ (35.0  *hh) 

D  =  (13.0  *  Ec  *  (1th  **  2))  /  (420.0  *  hh) 

E  =  (11.0  *  Ec  *  (1th  **  2))/ (210.0  *  hh) 

F  =  (9.0  *Ec*  1th)/ (70.0  *hh) 

G  =  (Ec  *  (1th  **  3))  /  (105.0  *  hh) 

H  =  (Ec*(lth**3))/(140.0*hh) 

Gc  =  Ec/eval 
Sil=Gc/hh 

Si2=Ec*hh/12.0/(eval**2) 


*  Input  Stiffiiess  Matrix 

Call  STIFF  1(A,B,K,L,P,R,S,T,W,Y,ST1) 
Call  STIFF2(C,D,E,F,G,H,ST2) 

Call  STIFF3(Sil,hh,hl,h2,lth,ST3) 

Call  STIFF4(Si2,hh,hl,h2,lth,ST4) 

DO  100,ROW=1,12 
DO  99,  COL=l,12 
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STIFFNESS(ROW,COL)=STl(ROW,COL)+ST2(ROW,COL)+ST3(ROW,COL) 

STIFFNESS(R0W,C0L)=STIFFNESS(R0W,C0L)+ST4(R0W,C0L) 

99  continue 

100  continue 

*  Add  Elements  FOR  K  MATRIX  and  M  MATRIX 

Mdim=6*(nel+1) 

DO  1000,  ROW=l,Mdim 
DO  1001,COL=l,Mdim 
Kmatrix(ROW,COL)  =  0 
Mmatrix(ROW,COL)  =  0 
1001  continue 
1000  CONTINUE 

*USE  Tins  COMMAND  TO  PRINT  THE  STIFFNESS  MATRIX  AND  ITS 
COMPONENTS 
If  (bit  .eq.  1)  then 

*  CALLWRRRN('K1 -uandw',  12, 12,Stl,  12,0) 

*  CALL  WRCRN  ('K2  -  phi',  12, 12,  St2, 12,  0) 

*  CALL  WRCRN  ('K3  -  si',  12, 12,  St3, 12, 0) 

*  CALL  WRCRN  ('K4  -  si  x',12,12,St4,12,0) 

*  Call  WRCRN('STIFF',12,12,STIFFNESS,12,0) 
end  if 


DO  1100,  IEL=1,NEL 
IES=  6*(IEL-1) 

DO  1102,  ROW  =1,12 
DO  1 101,  COL  =  1,12 
IS  =  ROW  +  lES 
JS  =  COL  +  lES 

Kmatrix(IS,JS)  =  Kmatrix(IS,JS)+STIFFNESS(ROW,COL) 
Mmatrix(IS,JS)  =  Mmatrix(IS,JS)+Mass(ROW,COL) 

1101  CONTINUE 

1102  CONTINUE 
1100  CONTINUE 


*  USE  THIS  COMMAND  TO  PRINT  THE  K  AND  M  MATRICES 

IF  (BIT  .EQ.  1)  THEN 
jed  =  6*(neH-l) 

*  Call  WRCRN('Kmatrix',jedjed,Kmatrix, 606,0) 

*  Call  WRRRN('Mmatrix',jed,jed,Mmatrix,606,0) 
end  if 


OMEGA=  FREQ'*2.0*3. 141 592654 
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*  Combine  Matrices 

DO  1200,  ROW=l,Mdim 
DO  1200,  COL=l,Mdim 

COMBINED(ROW,COL)=Kmatrix(ROW,COL)- 

((OMEGA**2)*Mmatrix(ROW,COL)) 

1200  CONTINUE 


♦Boundary  Conditions 

DO  1300,  ICR=l,NCONSTR 
NR=ICONNO(ICR) 

DO  1400,  ROW=l,Mdim 
COMBINED(ROW,NR)=0.0 
COMBINED(NR,ROW)=0.0 
1400  CONTINUE 

COMBINED(NR,NR)=l  .0 
1300  CONTINUE 

♦  Invert  Combined  Matrix 

CALL  LINCG  (Mdim,COMBINED,606,CINV,606) 


*  Multiply  F  and  CINV 


DO  2000,  ROW  =  l,Mdim 
Amatrix(ROW)  =  CINV(ROW,NMS)*Fo 
2000  CONTINUE 

Amatrix(NMS)=(  Amatrix(NMS)+Amatrix(NMS+3))/2 . 0 
Aconj=CONJG(Amatrix(NMS)) 

*  Calculate  Impedance 

Imped=1.0/OMEGA/SQRT(Amatrix(NMS)*Aconj) 

*  Print  Displacement  Vector  (For  given  values  of  fireq) 

*  If((FREQ.EQ.100.0).OR.(FREQ.EQ.160.0))THEN 

*  Print  *, ' ' 

*  Print  *, '  Peak  # ' 

*  Print  *,FREQ,Imped 

*  DO  200 1 ,  Nog=  1  ,Mdim 

*  Print  *,  Amatrix(Nog) 

*2001  continue 

*  end  if 

*  PRINT  IMPEDANCE  AND  FREQ 

Print  *,FREQ,Imped 
3000  CONTINUE 


END 


SUBROUTINE  STEEL(R,E) 

R=.000734 
E=  29000000.0 
END 

SUBROUTINE  OTHER(R,E) 

PRINT*, 'INPUT  DENSITY  (LB*SEC(squared)/INCH(to  the  fourth):' 
READ*,R 

PRINT*, 'INPUT  YOUNGS  MODULUS:' 

READ*,E 

END 

SUBROUTINE  NEOPRENE(0,R,EC,V) 

COMPLEX  EC 
*Poisson's  Ratio 
V-.49 

*  e 

E  =  2.0*(1.0  +  V) 

E  =  E 
*G 

Shearl  =  (2.71 828 18**(.026  *  LOG  (O)  +  4.754)) 

*  5 

Shear2  =  (2.7182818**(.17  *  LOG  (O)  -  2.705)) 

*  Modulus  for  Core 

EC  =  ((Shearl  *E),  (Shearl  *Shear2*E)) 

*  mass  density  of  core 

R=. 0001 15 
end 

SUBROUTINE  OCORE(R,EC) 

COMPLEX  EC 

PRINT  *,'INPUT  YOUNGS  MODULUS  (real,  imaginary):' 

READ  *,EC 

PRINT  *,'INPUT  DENSITY:' 

READ*,R 

end 

SUBROUTINE  STIFFl  (A,B,K,L,P,R,S,T,W,Y,STIFFNESS) 
INTEGER  ROW,  COL 
REAL  STIFFNESS(12,12),L,K 


*  Input  Stifihess  Matrix 
STIFFNESS(1, 1)  =  A 


STIFFNESS(1, 7)  = -A 
STIFFNESS(2, 2)  =  K 
STIFFNESS(2,  3)  =L 
STIFFNESS(2,  8)  =-K 
STIFFNESS(2, 9)=  L 
STIFFNESS(3, 3)  =  P 
STIFFNESS(3,  8)  =  -L 
STIFFNESS(3,  9)  =  R 
STIFFNESS(4, 4)  =B 
STIFFNESS(4,10)  =-B 
STIFFNESS(5,5)  =  S 
STIFFNESS(5,6)  =  T 
STIFFNESS(5,11)=-S 
STIFFNESS(5,12)=T 
STIFFNESS(6,  6)  =W 
STIFFNESS(6, 11)=-T 
STIFFNESS(6, 12)  =Y 
STIFFNESS(7,7>=A 
STIFFNESS(8,8)  =K 
STIFFNESS(8,9)  =-L 
STIFFNESS(9,9)  =P 
STIFFNESS(10,10)=B 
STIFFNESS(11,11)=S 
STIFFNESS(11,12)=-T 
STIFFNESS(12,12)=W 


DO  101,ROW=2,12 
IMl=ROW-l 
DO  101,COL=1,IM1 

STIFFNESS(ROW,COL)=STIFFNESS(COL,ROW) 
101  continue 
end 

SUBROUTINE  STIFF2(C,D,E,F,G,H, STIFFNESS) 
INTEGER  ROW,  COL 
COMPLEX  STIFFNESS(12, 12),C,D,E,F,G,H 
*  Input  Stiffness  Matrix 
STIFFNESS(2,  2)  =  C 
STIFFNESS(2, 3)  =E 
STIFFNESS(2,  5)  =  -C 
STIFFNESS(2, 6)  =-E 
STIFFNESS(2,  8)  =  F 
STIFFNESS(2, 9)  =-D 
STIFFNESS(2, 11)=-F 


STIFFNESS(2, 12)=  D 
STIFFNESS(3, 3)  =  G 
STIFFNESS(3,  5)  =-E 
STIFFNESS(3, 6)  =-G 
STIFFNESS(3,  8)  =D 
STIFFNESS(3, 9)  =  -H 
STIFFNESS(3, 11)=-D 
STIFFNESS(3, 12)  =H 
STIFFNESS(5,  5)  =C 
STIFFNESS(5,6)  =E 
STIFFNESS(5,8)  =-F 
STIFFNESS(5,9)  =  D 
STIFFNESS(5,11)=F 
STIFFNESS(5,12)=-D 
STIFFNESS(6, 6)  =G 
STIFFNESS(6,  8)  =  -D 
STIFFNESS(6, 9)  =H 
STIFFNESS(6, 11)  =D 
STIFFNESS(6, 12)  =-H 
STIFFNESS(8,8)=C 
STIFFNESS(8,9)  =-E 
STIFFNESS(8,11)=-C 
STIFFNESS(8,12)  =E 
STIFFNESS(9,9)  =G 
STIFFNESS(9,11)=E 
STIFFNESS(9,12)  =-G 
STIFFNESS(11, 11)=C 
STIFFNESS(11,  12)=-E 
STIFFNESS(12,12)  =G 


DO  102,  ROW=2,12 
IMl=ROW-l 
DO  102,  COL=  1,IM1 

STIFFNESS(ROW,COL)=STIFFNESS(COL,ROW) 
102  continue 
end 

SUBROUTINE  STIFF3(Si  1  ,hh,hl  ,h2,lth,STIFFNESS) 
INTEGER  ROW,  COL 
REAL  1th 

COMPLEX  STIFFNESS(12,12),Sil 
Si2=0.0 

*  Input  Stiffness  Matrix 

STIFFNESS(1, 1)  =  Sil*lth/3.0+Si2*1.0/lth 


STIFFNESS(1, 2)  =  -Sil*(hl+hh)/4.0 

STIFFNESS(1,  3)  =  Sil*(hl+hh)/24.0*lth+Si2*(hl+hh)/2.0/lth 

STIFFNESS(1, 4)  =  -Sil*lth/3.0-Si2*1.0/lth 

STIFFNESS(1,  5)  =  -Sil*(h2+hh)/4.0 

STIFFNESS(1, 6)  =Sil*(h2+hh)/24.0*lth+Si2*(h2+hh)/2.0/lth 

STIFFNESS(1, 7)  =  Sil*lth/6.0-Si2*1.0/lth 

STIFFNESS(1,  8)  =  Sil*(hl+hh)/4.0 

STIFFNESS(1, 9)  =  -Sil*(hl+hh)/24.0*lth-Si2*(hl+hh)/2.0/lth 

STIFFNESS(1, 10)  =  -Sil*lth/6.0+Si2*1.0/lth 

STIFFNESS(1, 11)  =  Sil*(h2+hh)/4.0 

STIFFNESS(1, 12)  =-Sil*(h2+hli)/24.0*lth-Si2*(h2+hh)/2.0/lth 

STIFFNESS(2,2)=((hl+hh)**2)*(Sil*6.0/5.0/lth+Si2*12.0/(lth**3)) 

STIFFNESS(2,2)=STIFFNESS(2,2)/4.0 

STIFFNESS(2,3)=((hl+hh)**2)/4.0*(Sil*1.0/10.0+Si2*6.0/(lth**2)) 
STIFFNESS(2, 4)  =  Sil*(hl+hh)/4.0 

STIFFNESS(2,5)=(hl+hh)*(h2+hh)*(Sil  *6.0/5.0/lth+Si2*  12.0/(lth**3)) 
STIFFNESS(2,5)=STIFFNESS(2,5)/4.0 

STIFFNESS(2,6)=(hl+hh)*(h2+hh)/4.0*(Sil*1.0/10.0+Si2*6.0/(lth**2)) 

STIFFNESS(2,6)=STIFFNESS(2,6) 

STIFFNESS(2, 7)  =  -Sil*(hl+hh)/4.0 

STIFFNESS(2,8)=((hl+hh)**2)*(-Sil  *6.0/5.0/lth-Si2*  12.0/(lth**3)) 
STIFFNESS(2,8)=STIFFNESS(2,8)/4.0 

STIFFNESS(2,9)=((hl+hh)**2)/4.0*(Sil*1.0/10.0+Si2*6.0/(lth**2)) 
STIFFNESS(2, 10)  =  Sil*(hl+hh)/4.0 

STIFFNESS(2,ll)=(h2+hh)/4.0*(-Sil*6.0/5.0/lth-Si2*12.0/(lth**3)) 
STIFFNESS(2,1 1)=STIFFNESS(2,1  l)*(hl+hh) 
STIFFNESS(2,12)=(hl+hh)*(h2+hh)*(Sil  *  1 .0/10.0+Si2*6.0/(lth**2)) 
STIFFNESS(2, 1 2)=STIFFNESS(2, 1 2)/4.0 
STIFFNESS(3,3)=((hl+hh)**2)/4.0*(Sil*2.0/15.0*lth+Si2*4.0/lth) 
STIFFNESS(3, 4)=-Sil  *(hl+hh)/24.0*lth-Si2*(hl+hh)/2.0/lth 
STIFFNESS(3,5)=STIFFNESS(2,12) 

STIFFNESS(3,6)=(hl+hh)*(h2+hh)/4.0*(Sil*2.0/15.0*lth+Si2*4.0/lth) 
STIFFNESS(3,6)=STIFFNESS(3 ,6) 

STIFFNESS(3,7)  =  -Sil*(hl+hh)/24.0*lth-Si2*(hl+hh)/2.0/lth 
STIFFNESS(3,8)=((hl+hh)**2)/4.0*(-Sil  *  1 .0/10.0-Si2*6.0/(lth**2)) 
STIFFNESS(3, 9)  =((hl+hh)**2)/4.0*(-Sil*lth/30.0+Si2*2.0/lth) 
STIFFNESS(3, 10)  =Sil*(hl4-hh)/24.0*lth+Si2*(hl+hh)/2.0/lth 
STIFFNESS(3,1  l)=(hl+hh)*(h2+hh)*(-Sil  *  1.0/10.0-Si2*6.0/(lth**2)) 
STIFFNESS(3,1 1)=STIFFNESS(3,1 1)/4.0 
STIFFNESS(3,12)=(hl+hh)*(h2+hh)/4.0*(-Sil*lth/30.0+Si2*2.0/lth) 
STIFFNESS(4, 4)  =Sil*lth/3.0+Si2*1.0/lth 
STIFFNESS(4,  5)  =  Sil*(h2+hh)/4.0 
STIFFNESS(4, 6)  =  -Si2*(li2+lih)/2.0/lth-Sil*(li2+hh)/24.0*lth 
STIFFNESS(4, 7)  =  -Sil*lth/6.0+Si2*1.0/lth 
STIFFNESS(4,  8)  =  -Sil*(hl+hh)/4.0 


STIFFNESS(4,  9)  =  Si2*(hl+hh)/2.0/lth+Sil*(hl+hh)/24.0*lth 
STIFFNESS(4, 10)  =+Sil*lth/6.0-Si2*1.0/lth 
STIFFNESS(4,  ll)=-Sil*(h2+hh)/4.0 
STIFFNESS(4, 12)  =  -STIFFNESS(4,6) 

STIFFNESS(5,5)=((h2+hh)**2)*(Sil  *6.0/5.0/lth+Si2*  12.0/(lth**3)) 
STIFFNESS(5,5)=STIFFNESS(5,5)/4.0 

STIFFNESS(5,6)=((h2+hh)**2)/4.0*(Si  1  *  1 .0/1 0.0+Si2*6.0/(lth**2)) 
STIFFNESS(5, 7)  =-Sil*(h2+hh)/4.0 
STIFFNESS(5,  8)  =  STIFFNESS(2,1 1) 

STIFFNESS(5,  9)  =  STIFFNESS(2,12) 

STIFFNESS(5, 10)  =Sil*(h2+hh)/4.0 

STIFFNESS(5,1  l)=((h2+hh)**2)*(-Sil  *6.0/5.0/lth-Si2*  12.0/(lth**3)) 
STIFFNESS(5,1 1)=  STIFFNESS(5,1 1)/4.0 
STIFFNESS(5,12)=((h2+hh)**2)/4.0*(Sil*1.0/10.0+Si2*6.0/(lth**2)) 
STIFFNESS(6,6)=((h2+hh)**2)/4.0*(Sil  *2.0/1 5.0*lth+Si2*4.0/lth) 
STIFFNESS(6,7)=-Si2*(h2+hh)/2.0/lth-Sil*(h2+hh)/24.0*lth 
STIFFNESS(6,  8)  =  STIFFNESS(3,1 1) 

STIFFNESS(6,  9)  =  STIFFNESS(3,12) 

ST1FFNESS(6, 10)  =  Si2*(h2+hh)/2.0/lth+Sil*(h2+hli)/24.0*lth 
STIFFNESS(6, 1 1  )=((h2+hh)*  *2)/4.0*(-Si  1  *  1 .0/1 0.0-Si2*6.0/(lth*  *2)) 
STIFFNESS(6,1 1)=STIFFNESS(6,1 1) 

STIFFNESS(6, 12)  =((h2+hh)**2)/4.0*(-Sil*lth/30.0+Si2*2.0/lth) 

STIFFNESS(7, 7)  =  Sil*lth/3.0+Si2*1.0/lth 

STIFFNESS(7,  8)  =  Sil*(hl+hh)/4.0 

STIFFNESS(7, 9)  =  Si2*(hl+hh)/2.0/lth+Sil*(hl+hh)/24.0*lth 

STIFFNESS(7, 10)  =  -Sil*lth/3.0-Si2*1.0/lth 

STIFFNESS(7, 1 1)  =  Sil  *(h2+hh)/4.0 

STIFFNESS(7, 12)  =  Si2*(h2+hh)/2.0/lth+Sil*(h2+hh)/24.0*lth 

STIFFNESS(8,  8)  =  STIFFNESS(2,2) 

STIFFNESS(8,  9)  =  STIFFNESS (3 ,8) 

STIFFNESS(8, 10)  =-Sil*(hl+hh)/4.0 
STIFFNESS(8, 1 1)  =  STIFFNESS(2,5) 

STIFFNESS(8, 12)  =STIFFNESS(3,11) 

STIFFNESS(9, 9)  =STIFFNESS(3,3) 

STIFFNESS(9, 10)=-Si2*(hl+hh)/2.0/lth-Sil*(hl+hh)/24.0*lth 
STIFFNESS(9, 11)  =STIFFNESS(8,12) 

STIFFNESS(9, 12)  =STIFFNESS(3,6) 

STIFFNESS(10, 10)  -Sil*lth/3.0+Si2*1.0/lth 
STIFFNESS(10, 1 1)  =  -Sil*(h2+hh)/4.0 
STIFFNESS(10, 12)  =-Si2*(h2+hh)/2.0/lth-Sil*(h2+hh)/24.0*lth 
STIFFNESS(11, 11)=STIFFNESS(5,5) 

STIFFNESS(1 1, 12)  =STIFFNESS(6,1 1) 

STIFFNESS(12, 12)  =  STIFFNESS(6,6) 


DO  103,  ROW=2,12 


IM1=R0W-1 
DO  103,COL=1,IM1 

STIFFNESS(ROW,COL)=STIFFNESS(COL,ROW) 

1 03  continue 
end 

SUBROUTINE  STIFF4(Si2,hh,hl  ,h2, 1th, STIFFNESS) 

INTEGER  ROW,  COL 
REAL  1th 

COMPLEX  STIFFNESS(12,12),Si2 
Si  1=0.0 

*  Input  Stiffness  Matrix 

STIFFNESS(1, 1)  =  Sil*lth/3.0+Si2*1.0/lth 

STIFFNESS(1, 2)  =  -Sil*(hl+hh)/4.0 

STIFFNESS(1, 3)  =  Sil  *(hl+hh)/24.0*lth+Si2*(hl+hh)/2.0/lth 

STIFFNESS(1, 4)  =  -Sil*lth/3.0-Si2*1.0/lth 

STIFFNESS(1,  5)  =  -Sil*(h2+hh)/4.0 

STIFFNESS(1, 6)  =Sil*(h2+hh)/24.0*lth+Si2*(h2+hh)/2.0/lth 

STIFFNESS(1, 7)  =  Sil*lth/6.0-Si2*1.0/lth 

STIFFNESS(1, 8)  =  Sil*(hl+hh)/4.0 

STIFFNESS(1, 9)  =  -Sil*(hl+hh)/24.0*lth-Si2*(hl+hh)/2.0/lth 

STIFFNESS(1, 10)  =  -Sil*lth/6.0+Si2*1.0/lth 

STIFFNESS(1,  ll)  =  Sil*(h2+hh)/4.0 

STIFFNESS(1, 12)  =-Sil*(h2+hh)/24.0*lth-Si2*(h2+hh)/2.0/lth 

STIFFNESS(2,2)=((hl+hh)**2)*(Sil*6.0/5.0/lth+Si2*12.0/(lth**3)) 

STIFFNESS(2,2)=STIFFNESS(2,2)/4.0 

STIFFNESS(2,3)=((hl+hh)**2)/4.0*(Sil*1.0/10.0+Si2*6.0/(lth**2)) 
STIFFNESS(2, 4)  =  Sil*(hl+hh)/4.0 

STIFFNESS(2,5)=(hl+hh)*(h2+hh)*(Sil*6.0/5.0/lth+Si2*12.0/(lth**3)) 

STIFFNESS(2,5)=STIFFNESS(2,5)/4.0 

STIFFNESS(2,6)=(hl  +hh)*(h2+hh)/4.0*(Si  1  ♦  1 .0/1 0.0+Si2*6.0/(lth*  *2)) 
STIFFNESS(2,6)=STIFFNESS(2,6) 

STIFFNESS(2,  7)  =  -Sil*(hl+hh)/4.0 

STIFFNESS(2,8)=((hl+hh)**2)*(-Sil  *6.0/5. 0/lth-Si2*  12.0/(lth**3)) 
STIFFNESS(2,8)=STIFFNESS(2,8)/4.0 

STIFFNESS(2,9)=((hl+hh)**2)/4.0*(Sil  *  1 .0/10.0+Si2*6.0/(lth**2)) 
STIFFNESS(2, 10)  =  Sil*(hl+hh)/4.0 

STIFFNESS(2,ll)=(h2+hh)/4.0*(-Sil*6.0/5.0/lth-Si2*12.0/(lth**3)) 

STIFFNESS(2, 1 1  )=STIFFNESS(2, 1 1  )*(hl  +hh) 

STIFFNESS(2,12)=(hl+hh)*(h2+hh)*(Sil*1.0/10.0+Si2*6.0/(lth**2)) 

STIFFNESS(2,12)=STIFFNESS(2,12)/4.0 

STIFFNESS(3,3)=((hl+hh)**2)/4.0*(Sil*2.0/15.0*lth+Si2*4.0/lth) 

STIFFNESS(3, 4)=-Sil  *(hl+hh)/24.0*lth-Si2*(hl+hh)/2.0/lth 

STIFFNESS(3,5)=STIFFNESS(2,12) 

STIFFNESS(3,6)=(hl+hh)*(h2+hh)/4.0*(Sil*2.0/15.0*lth+Si2*4.0/lth) 


STIFFNESS(3,6)=STIFFNESS(3,6) 

STIFFNESS(3,7)  =  -Sil*(hl+hh)/24.0*lth-Si2*(hl+hh)/2.0/lth 
STIFFNESS(3,8)=((hl+hh)**2)/4.0*(-Sil  *  1 .0/10.0-Si2*6.0/(lth**2)) 
STIFFNESS(3, 9)  =((hl+hh)**2)/4.0*(-Sil*lth/30.0+Si2*2.0/lth) 
STIFFNESS(3, 10)  =Sil*(hl+hh)/24.0*lth+Si2*(hl+hh)/2.0/lth 
STIFFNESS(3,1  l)=(hl+hh)*(h2+hh)*(-Sil  *  1.0/10.0-Si2*6.0/(lth**2)) 
STIFFNESS(3,1 1)=STIFFNESS(3,1 1)/4.0 

STIFFNESS(3,12)=(hl+hh)*(h2+hh)/4.0*(-Sil*lth/30.0+Si2*2.0/lth) 
STIFFNESS(4, 4)  =Sil*lth/3.0+Si2*1.0/Ith 
STIFFNESS(4,  5)  =  Sil*(h2+hh)/4.0 

STIFFNESS(4, 6)  =  -Si2*(h2+hh)/2.0/lth-Sil*(h2+hh)/24.0*lth 

STIFFNESS(4, 7)  =  -Sil*lth/6.0+Si2*1.0/lth 

STIFFNESS(4,  8)  =  -Sil*(hl+hh)/4.0 

STIFFNESS(4, 9)  =  Si2*(hl+hh)/2.0/lth+Sil*(hl+hh)/24.0*lth 

STIFFNESS(4, 10)  =+Sil*lth/6.0-Si2*1.0/lth 

STIFFNESS(4,  ll)=-Sil*(h2+hh)/4.0 

STIFFNESS(4, 12)  =  -STIFFNESS(4,6) 

STIFFNESS(5,5)=((h2+hh)**2)*(Sil*6.0/5.0/lth+Si2*12.0/(lth**3)) 

STIFFNESS(5,5)=STIFFNESS(5,5)/4.0 

STIFFNESS(5,6)=((h2+hh)**2)/4.0*(Sil*1.0/10.0+Si2*6.0/(lth**2)) 
STIFFNESS(5, 7)  =-Sil*(h2+hh)/4.0 
STIFFNESS(5,  8)  =  STIFFNESS(2,1 1) 

STIFFNESS(5,  9)  =  STIFFNESS(2,12) 

STIFFNESS(5, 10)  =Sil*(h2+hh)/4.0 

STIFFNESS(5, 1 1  )=((h2+hh)*  *2)*(-Sil  *6.0/5.0/lth-Si2*  1 2.0/(lth*  *3)) 
STIFFNESS(5,1 1)=  STIFFNESS(5,1 1)/4.0 

STIFFNESS(5,12)=((h2+hh)**2)/4.0*(Sil*1.0/10.0+Si2*6.0/(lth**2)) 
STIFFNESS(6,6)=((h2+hh)**2)/4.0*(Sil*2.0/15.0*lth+Si2*4.0/lth) 
STIFFNESS(6,7)=-Si2*(h2+hh)/2.0/lth-Sil*(h2+hh)/24.0*lth 
STIFFNESS(6,  8)  =  STIFFNESS(3,1 1) 

STIFFNESS(6, 9)  =  STIFFNESS(3,12) 

STIFFNESS(6, 10)  =  Si2*(h2+hh)/2.0/lth+Sil*(h2+hh)/24.0*lth 
STIFFNESS(6, 1 1  )=((h2+hh)*  *2)/4.0*  (-Si  1  *  1 .0/1 0.0-Si2*  6.0/(lth*  *2)) 
STIFFNESS(6,1 1)=STIFFNESS(6,1 1) 

STIFFNESS(6, 12)  =((h2+hh)**2)/4.0*(-Sil*lth/30.0+Si2*2.0/lth) 

STIFFNESS(7, 7)  =  Sil*lth/3.0+Si2*1.0/lth 

STIFFNESS(7,  8)  =  Sil*(hl+hh)/4.0 

STIFFNESS(7, 9)  =  Si2*(hl+hh)/2.0/lth+Sil*(hl+hh)/24.0*lth 

STIFFNESS(7, 10)  =  -Sil*lth/3.0-Si2*1.0/lth 

STIFFNESS(7, 11)  =  Sil*(h2+hh)/4.0 

STIFFNESS(7, 12)  =  Si2*(h2+hh)/2.0/lth+Sil*(h2+hh)/24.0*lth 

STIFFNESS(8,  8)  =  STIFFNESS(2,2) 

STIFFNESS(8,  9)  =  STIFFNESS (3, 8) 

STIFFNESS(8, 10)  =-Sil*(hl+hh)/4.0 
STIFFNESS(8, 1 1)  =  STIFFNESS(2,5) 


STIFFNESS(8, 12)  =STIFFNESS(3,1 1) 

STIFFNESS(9, 9)  =STIFFNESS(3,3) 

STIFFNESS(9, 10)=-Si2*(hl+hh)/2.0/lth-Sil*(hl+hh)/24.0*lth 
STIFFNESS(9, 11)  =STIFFNESS(8,12) 

STIFFNESS(9, 12)  =STIFFNESS(3,6) 

STIFFNESS(10, 10)  =Sil*lth/3.0+Si2*1.0/lth 
STIFFNESS(10, 1 1)  =  -Sil*(h2+hh)/4.0 
STIFFNESS(10, 12)  =-Si2*(h2+hh)/2.0/lth-Sil*(h2+hh)/24.0*lth 
STIFFNESS(1 1, 1 1)  =STIFFNESS(5,5) 

STIFFNESS(11, 12)  =STIFFNESS(6,11) 

STIFFNESS(12, 12)  =  STIFFNESS(6,6) 

DO  104,  ROW=2,12 
IMl=ROW-l 
DO  104,  COL=  1,IM1 

STIFFNESS(ROW,COL)=STIFFNESS(COL,ROW) 

104  continue 
end 


SUBROUTINE  Massl  (hh,hl  ,h2,lth,rr,eval,Mass) 

INTEGER  ROW,  COL 
REAL  1th,  Mass(12,12) 

U 1 6=rr*(6.0*hh*(hl +hh)-5.0*eval*(lth*  *2)) 
U16=U16*((hh/eval)**2)/1440.0/lth 

U55=42.0*hh*hl*((hh**2)-eval*(lth**2)+hh*hl/2.0)+2L0*(hh**4) 

U55-U55+104.0*(eval**2)*(lth**4)-42.0*eval*((hh*lth)**2) 

U55=U55*hh*rr/840.0/(lth**3)/(eval**2) 

U56=126.0*hh*hl*((hh**2)-eval*(lth**2)+hh*hl/2.0)+63.0*(hh**4) 

U56=U56+88.0*(eval**2)*(lth**4)-126.0*eval*((hh*lth)**2) 

U56=U56*hh*rr/5040.0/(lth**2)/(eval**2) 

U57=14.0*hh*hl*((hh**2)-eval*(lth**2)+hh*hl/2.0)+7.0*(hh**4) 

U57=U57-12.0*(eval**2)*(lth**4)-14.0*eval*((hh*lth)**2) 

U57=U57*(-hh)*rr/280.0/(lth**3)/(eval**2) 

U58=126.0*hh*hl*((hh**2)-eval*(lth**2)/6.0+hh*hl/2.0)+63.0*(hh**4) 

U58=U58-52.0*(eval**2)*(lth**4)-2L0*eval*((hh*lth)**2) 

U58=U58*hh*rr/5040.0/(lth**2)/(eval**2) 

U66=42.0*hh*hl*((hh**2)-eval*(Ith**2)/3.0+hh*hl/2.0)+21.0*(hh**4) 

U66=U66+8.0*(eval**2)*(lth**4)-14.0*eval*((hh*lth)**2) 

U66=U66*hh*rr/2520.0/(lth)/(eval**2) 

U00=42.0*hh*h2*((hh**2)-eval*(lth**2)/3.0+hh*h2/2.0)+21.0*(hh**4) 

U00=U00+8.0*(eval**2)*(lth**4)-14.0*eval*((hh*lth)**2) 

U00=U00*hh*rr/2520.0/(lth)/(eval**2) 

U10=126.0*hh*h2*(eval*(lth**2)/6.0-(hh**2)-hh*h2/2.0)-63.0*(hh**4) 

U10=U10+52.0*(eval**2)*(lth**4)+21.0*eval*((hh*lth)**2) 
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U1 0=U1 0*hh*rr/5040.0/(lth*  *2)/(eval*  *2) 

U68=42.0*hh*hl*((hh**2)+eval*(lth**2)/6.0+hh*hl/2.0)+21.0*(hh**4) 

U68=U68-12.0*(eval**2)*(lth**4)+7.0*eval*((hh*lth)**2) 

U68=U68*hh*rr/5040.0/(lth)/(eval**2) 

U99=42.0*hh*h2*((hh**2)-eval*(lth**2)+hh*h2/2.0)+21.0*(hh**4) 

U99=U99+104.0*(eval**2)*(lth**4)-42.0*eval*((hh*lth)**2) 

U99=U99*hh*rr/840.0/(lth**3)/(eval**2) 

U90=126.0*hh*h2*((hh**2)-eval*(lth**2)+hh*h2/2.0)+63.0*(hh**4) 
U90=U90+88.0*(eval**2)*(lth**4)-126.0*eval*((hh*lth)**2) 
U90=U90*hh*ir/5040.0/(lth*  *2)/(eval*  *2) 
U91=14.0*hh*h2*((hh**2)-eval*(lth**2)+hli*h2/2.0)+7.0*(hh**4) 
U91=U91-12.0*(eval**2)*(lth**4)-14.0*eval*((hh*lth)**2) 

U9 1=U9 1  *(-hh)*rr/280.0/(lth**3)/(eval*  *2) 

U1  l=42.0*hh*h2*((hh**2)+eval*(lth**2)/6.0+hh*h2/2.0)+21 .0*(hh**4) 

U1 1=U1  l-12.0*(eval**2)*(lth**4)+7.0*eval*((hh*lth)**2) 

U1 1=U1  l*hh*iT/5040.0/(lth)/(eval**2) 

U59=21 .0*hh*(hh*hl  *h2-eval*hl  *(lth*  *2)+(hh*  *3)+h2*hh*hh+hl  *hh*hh) 

U59=U59+52.0*(eval**2)*(lth**4)-21.0*hh*eval*h2*(lth**2) 

U59=(U59-42.0*eval*((hh*lth)**2))*hh*rr/840.0/(Ith**3)/(eval**2) 

U69=126.0*hh*(hh*hl*h2-eval*h2*(lth**2)/6.0+(hh**3)+h2*hh*hh) 

U69=U69+88.0*(eval**2)*(lth**4)-231.0*hh*eval*hl*(lth**2) 

U69=U69+126.0*(hh**3)*hl-252.0*eval*((hh*lth)**2) 

U69=U69*hh*rr/l  0080.0/(lth**2)/(eval*  *2) 

U79=7.0*hh*(eval*hl*(lth**2)-hh*hl*h2-(hh**3)-h2*hh*hh-hl*hh*hh) 

U79=U79+6.0*(eval**2)*(lth**4)+7.0*hh*eval*h2*(lth**2) 

U79=(U79+14.0*eval*((hh*lth)**2))*hh*rr/280.0/(lth**3)/(eval**2) 

U89=l  26.0*hh*(hh*hl  *h2-eval*h2*(lth*  *2)/6.0+(hh*  *3)+h2*hh*hh) 

U89-U89-52.0*(eval**2)*(lth**4)-21.0*hh*eval*hl*(lth**2) 

U89=U89+126.0*(hh**3)*hl-42.0*eval*((hh*lth)**2) 

U89=U89*hh*rr/l  0080.0/(lth**2)/(eval*  *2) 

U50=126.0*hh*(hh*hl  *h2-eval*hl  *(lth**2)/6.0+(hh**3)+h2*hh*hh) 
U50=U50+88.0*(eval**2)*(lth**4)-231.0*hh*eval*h2*(lth**2) 
U50=U50-252.0*eval*((hh*lth)**2)+126.0*(hh**3)*hl 
U50=U50*hh*rr/10080.0/(lth**2)/(eval**2) 

U60=21 .0*hh*(bh*hl  *h2-eval*h2*(lth**2)/3.0+(hh**3)+h2*hh*hh) 
U60=U60+4.0*(eval**2)*(lth**4)-7.0*hh*eval*hl*(lth**2) 
U60=U60+21.0*(hh**3)*hl-14.0*eval*((hh*lth)**2) 
U60=U60*hh*rr/2520.0/(lth)/(eval**2) 

U80=42.0*hh*(eval*hl*(ith**2)/6.0+hh*hl*h2+(hh**3)+h2*hh*hh) 

U80=U80-12.0*(eval**2)*(lth**4)+7.0*hh*eval*h2*(lth**2) 

U80=U80+14.0*eval*((hh*lth)**2)+42.0*(hh**3)*hl 

U80=U80*hh*rr/10080.0/(lth)/(eval**2) 

Mass(l ,  1  )=rr*(hh*  *3)/l  20.0/lth/(eval**2) 
Mass(l,7)=-rr*(hh**3)/120.0/lth/(eval**2) 

Mass(l  ,4)=-rr*(hh**3)/l  20.0/lth/(eval*  *2) 


Mass(l,10)=rr*(hh**3)/120.0/lth/(eval**2) 
Mass(l  ,2)=-rr*(hh**2)/48.0/eval 
Mass(l,3)=U16 

Mass(l,8)=-rr*(hh**2)/48.0/eval 

Mass(l,9)=-U16 

Mass(  1 ,5)=-rr*(hh’'' *2)/48.0/eval 

Mass(l,6)=U16 

Mass(l,l  l)=-rr*(]ih**2)/48.0/eval 
Mass(l,12)=-U16 

Mass(7,7)=rr*(hh*  *3)/l  20.0/lth/(eval*  *2) 

Mass(4,7)=rr*(hh*  *3)/l  20.0/lth/(eval*  *2) 

Mass(7,l  0)=-iT*(hh**3)/120.0/lth/(eval*  *2) 

Mass(2,7)=rr*(hh**2)/48.0/eval 

Mass(3,7)=-U16 

Mass(7,8)=rT*(hh**2)/48.0/eval 

Mass(7,9)=U16 

Mass(5,7)=rr*(hh*  *2)/48.0/eval 
Mass(6,7)=-U16 

Mass(7,l  l)=rr*(hh**2)/48.0/eval 
Mass(7,12)=U16 

Mass(4,4)=rr*(hh**3)/120.0/lth/(eval**2) 
Mass(4, 1 0)=-rr*(hh**3)/l  20.0/lth/(eval*  *2) 
Mass(2,4)=iT*(hh*  *2)/48 .0/eval 
Mass(3,4)=-U16 
Mass(4,8)=rr*(hli**2)/48.0/eval 
Mass(4,9)=U16 

Mass(4,5)=Tr*(hh*’''2)/48.0/eval 

Mass(4,6)=-U16 

Mass(4,l  l)=rr*(hh**2)/48.0/eval 
Mass(4,12)=U16 

Mass(l  0, 1 0)=rr*(hh**3)/l  20.0/lth/(eval*  *2) 
Mass(2, 1 0)=-rr*(hh*  *2)/48.0/eval 
Mass(3,10)=U16 

Mass(8, 1 0)=-rr*(hh*  *2)/48.0/eval 

Mass(9,10)=-U16 

Mass(5,l  0)=-rr*(hh**2)/48.0/eval 

Mass(6,10)=U16 

Mass(10,l  l)=-rr*(hh**2)/48.0/eval 

Mass(10,12)=-U16 

Mass(2,2)=U55 

Mass(2,3)-U56 

Mass(2,8)=U57 

Mass(2,9)=U58 

Mass(2,5)=U59 

Mass(2,6)=U50 
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Mass(2,ll)=U79 
Mass(2,12)=U89 
Mass(3,3)=U66 
Mass(3,8)=-U58 
Mass(3,9)=U68 
Mass(3,5)=U69 
Mass(3,6)=U60 
Mass(3,ll)=-U89 
Mass(3,12)-U80 
Mass(8,8)=U55 
Mass(8,9)=-U56 
Mass(5,8)=U79 
Mass(6,8)=-U89 
Mass(8,ll)=U59 
Mass(8,12)=-U50 
Mass(9,9)=U66 
Mass(5,9)=U89 
Mass(6,9)=U80 
Mass(9,ll)=-U69 
Mass(9,12)=U60 
Mass(5,5)=U99 
Mass(5,6)=U90 
Mass(5,ll)=U91 
Mass(5,12)=-U10 
Mass(6,6)=U00 
Mass(6,ll)=U10 
Mass(6,12)=Ull 
Mass(ll,ll)=U99 
Mass(ll,12)=-U90 
Mass(12,12)=U00 
DO  50,  ROW=2,12 
IMl=ROW-l 
DO  50,  COL=  1,IM1 
Mass(row,col)=Mass(coI,row) 

50  continue 
end 

SUBROUTINE  Mass2(hh,rr,eval,hl  ,h2,lth,  Alpha,Mass) 

INTEGER  ROW,  COL 

REAL  Mas(12,12),Mass(12,12),lth 

Wll  =1.0/3 .0*hh*rr 

W21=-l  .0/96.0*(hh**3)/eval*rr 

W22  =1.0/768.0*(hh**5)/(eval**2)*rr 

W51=1.0/24.0*hh*(4.0*hl+hh)*rr 

W53  =1. 0/24.0 W(2.0*hl+hh)*rr 


W83  =  -1.0/24.0*hh*(4.0*h2+hh)*rr 

W81  =-1.0/24.0*hh*(2.0*h2+hh)*rr 

W52  =  -1.0/192.0*(hh**3)*hl/eval*rr 

W55-1.0/120.0*hh*(10.0*(hl**2)+5.0*hl*hh+(hh**2))*rr 

W82  =  -1.0/192.0*(hh**3)*h2/eval*rr 

W61=-1.0/192.0*(hl+hh)*(hh**3)*rr/eval 

W62=1.0/1536.0*(hh**5)*(hl+hh)/(eval**2)*rr 

W65=-l  .0/384.0*(hl+hh)*(hh**3)*hl  *rr/eval 

W71=1.0/960.0*(hH-hh)*(hh**3)*rr/eval 

W72=-1.0/7680.0*(hh**5)*(hl+hh)/(eval**2)*rr 

W75=1.0/1920.0*(hl+hh)*(hh**3)*hl*rr/eval 

W85=-1.0/240.0*hh*(2.0*(hh**2)+5.0*hh*(hl+h2)+10.0*hl*h2)*rr 

W66=1.0/3072.0*(hh**5)*((hl+hh)**2)/(eval**2)*rr 

W76=-1.0/15360.0*(hh**5)*((hl+hh)**2)/(eval**2)*rr 

W77=l  .0/645 1 2.0*(hh*  *5)*((hl +hh)*  *2)/(eval*  *2)*rr 

W86=-1.0/384.0*(hl+hh)*(hh**3)*h2*rr/evaI 

W87=l  .0/1 920.0*(hl+hh)*(hh*  *3)*h2*rr/eval 

W88=1.0/120.0*hh*((hh**2)+5.0Wh2+10.0*(h2**2))*rr 

W91=-1.0/192.0*(h2+hh)*(hh**3)*rr/eval 

W92=1.0/1536.0*(hh**5)*(h2+hh)/(eval**2)*rr 

W95=-l  .0/384.0*(h2+hh)*(hh**3)*hl  *rr/eval 

W01=1.0/960.0*(h2+hh)*(hh**3)*rr/eval 

W02=-1.0/7680.0*(hh**5)*(h2+hh)/(eval**2)*rr 

W05=l  .0/1 920.0*(h2+hh)*(hh*  *3)*hl  *rr/eval 

W96=1.0/3072.0*(hh**5)*(hl+hh)*(h2+hh)/(eval**2)*rr 

W97=-1.0/15360.0*(hh**5)*(hl+hh)*(h2+hh)/(eval**2)*rr 

W98=-1.0/384.0*(h2+hh)*(hh**3)*h2*rr/eval 

W99=1.0/3072.0*(hh**5)*((h2+hh)**2)/(eval**2)*rr 

W07=l  .0/645 1 2.0*(hh*  *5)*(hl +hh)*(h2+hh)/(eval*  *2)*rr 

W08=l  .0/1 920.0*(h2+hh)*(hh*  *3)*h2*rr/eval 

W09=-1.0/15360.0*(hh**5)*((h2+hh)**2)/(eval**2)*rr 

W00=1 .0/645 1 2.0*(hh*  *  5)*((h2+hh)**2)/(eval*  *2)*rr 

Mas(l ,  1  )=(- 1 .0/3 .0/lth*(-Wl  1  *(lth*  *2)-3 .0*  W22+3.0*  W21  *lth)) 

Mas(l  ,2)=1 .0/6.0/lth*(-6.0*  W22+W1 1  *(lth**2)) 
Mas(l,3)=1.0/6.0/lth*(Wll*(lth**2)+6.0*W21*lth-6.0*W22) 
Mas(l,4)=1.0/6.0/lth*(Wll/2.0*(lth**2)+W22*6.0) 
Mas(l,5)=ALPHA*lth*(-2.0*W72+lth*W71)-lth*W51+2.0*(W52-W61) 
Mas(l  ,5)=1 .0/2.0/lth*Mas(l  ,5) 

Mas(l  ,6)=1 .0/12.0/lth*(12.0*(W62-W6 1  *lth)+W5 1  *(lth**2)) 
Mas(l,7)=ALPHA*lth*(-2.0*W72+lth*W71)+lth*W51+2.0*(W61-W52) 
Mas(l  ,7)=1 .0/2.0/lth*Mas(l,7) 

Mas(l,8)=-1 .0/12.0/lth*(12.0*W62+W5 1  *(lth**2)) 
Mas(l,9)=ALPHA*lth*(-2.0*W02+lth*W01)-lth*W81+2.0*(W82-W91) 
Mas(l  ,9)=1 .0/2.0/lth*Mas(l  ,9) 

Mas(l,10)=1.0/12.0/lth*(12.0*(W92-W91*lth)+W81*(lth**2)) 


Mas(l,ll)=ALPHA*lth*(-2.0*W02+lth*W01)+lth*W81+2.0*(W91-W82) 
Mas(l,l  l)=1.0/2.0/lth*Mas(l,l  1) 
Mas(l,12)=-1.0/12.0/Ith*(12.0*W92+W81*(lth**2)) 
Mas(2,2)=(1.0/3.0/lth*(Wll*(lth**2)+3.0*W22+3.0*W21*lth)) 
Mas(2,3)=Mas(l,4) 

Mas(2,4)=l  .0/6.0/lth*(Wl  1  *(lth**2)-6.0*  W21  *lth-6.0*  W22) 
Mas(2,5)=ALPHA*lth*(2.0*W72+lth*W71)-lth*W51-2.0*(W52-W61) 
Mas(2,5)=l  .0/2.0/lth*Mas(2,5) 

Mas(2,6)=Mas(l,8) 

Mas(2,7)=ALPHA*lth*(2.0*  W72+lth*  W7 1  )+lth*  W5 1+2.0*(W52-W6 1 ) 
Mas(2,7)=l  .0/2.0/lth*Mas(2,7) 

Mas(2,8)=1.0/12.0/lth*(12.0*(W62+W61*lth)+W51*(lth**2)) 
Mas(2,9)=ALPHA*lth*(2.0*W02+lth*W01)-lth*W81+2.0*(W91-W82) 
Mas(2,9)=l  .0/2.0/lth*Mas(2,9) 

Mas(2,10)=Mas(l,12) 

Mas(2, 1 1  )=ALPHA*lth*(2.0*  W02+lth*  W01)+lth*  W8 1  -2.0*(W9 1 -W82) 
Mas(2,l  1)-1 .0/2.0/lth*Mas(2,l  1) 

Mas(2,12)=1.0/12.0/lth*(12.0*(W92+W91*lth)+W81*(lth**2)) 

Mas(3,3)=Mas(l,l) 

Mas(3,4)=Mas(l,2) 

Mas(3,5)=ALPHA*lth*(-2.0*W72+lth*W71)+lth*W53+2.0*(W52-W61) 
Mas(3,5)=-1 .0/2.0/lth*Mas(3,5) 

Mas(3,6)=1.0/12.0/lth*(12.0*(-W62+W61*lth)+W53*(lth**2)) 
Mas(3,7)=ALPHA*lth*(-2.0*W72+lth*W71)-lth*W53+2.0*(W61-W52) 
Mas(3,7)=-1 .0/2.0/lth*Mas(3,7) 
Mas(3,8)=1.0/12.0/lth*(12.0*W62-W53*(lth**2)) 
Mas(3,9)=ALPHA*lth*(-2.0*W02+lth*W01)+lth*W83+2.0*(W82-W91) 
Mas(3,9)=-1 .0/2.0/Ith*Mas(3,9) 

Mas(3,l  0)=- 1 .0/1 2.0/lth*(l  2.0*(W92-W91  *lth)-W83  *(lth*  *2)) 
Mas(3,ll)=ALPHA*lth*(-2.0*W02+lth*W01)-lth*W83+2.0*(W91-W82) 
Mas(3,l  l)=-1.0/2.0/lth*Mas(3,ll) 
Mas(3,12)=1.0/12.0/lth*(12.0*W92-W83*(lth**2)) 

Mas(4,4)=Mas(2,2) 

Mas(4,5)=ALPHA*lth*(2.0*W72+lth*W71)+lth*W53-2.0*(W52-W61) 
Mas(4,5)=-1 .0/2.0/lth*Mas(4,5) 

Mas(4,6)=Mas(3,8) 

Mas(4,7)=ALPHA*lth*(2.0*W72+lth*W71)-lth*W53+2.0*(W52-W61) 
Mas(4,7)=-1 .0/2.0/lth*Mas(4,7) 

Mas(4,8)=-1 .0/12.0/lth*(12.0*(W62+W61  *lth)-W53*(lth**2)) 
Mas(4,9)=ALPHA*lth*(2.0*W02+lth*W01)+lth*W83+2.0*(W91-W82) 
Mas(4,9)=-1 .0/2.0/lth*Mas(4,9) 

Mas(4, 1 0)=Mas(3 , 1 2) 

Mas(4,ll)=ALPHA*lth*(2.0*W02+lth*W01)-lth*W83-2.0*(W91-W82) 
Mas(4,l  1)=-1 .0/2.0/lth*Mas(4,l  1) 

Mas(4,12)=-1 .0/12.0/lth*(12.0*(W92+W91  *lth)-W83*(lth**2)) 
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Mas(5,5)=ALPHA*(lth**3)*5.0*(-2.0*W75+lth*W77*ALPHA)+60.0*W66 

Mas(5,5)=1.0/5.0/(lth**3)*(Mas(5,5)+6.0*W55*(lth**2)) 

Mas(5,6)=- 1 .0/1 0.0*(1 0.0*  W76*  ALPHA- W55-W66*60.0/(lth*  *2)) 

Mas(5,7)=(ALPHA**2)*W77*(lth**4)*5.0-W55*6.0*(lth**2)-W66*60.0 

Mas(5,7)=1.0/5.0/(lth**3)*Mas(5,7) 

Mas(5,8)=1.0/10.0*(10.0*W76*ALPHA+W55+W66*60.0/(lth**2)) 

Mas(5,10)=W95/lth-ALPHA*W97-W86/lth+W85/10.0+6.0*W96/(lth**2) 

Mas(5,9)=ALPHA*(ALPHA*lth*W07-W05-W87)+6.0*W85/lth/5.0 

Mas(5,9)=Mas(5,9)+60.0*W96/(lth**3)/5.0 

Mas(5, 1 1  )=ALPHA*(ALPHA*Ith*  W07-W05+W87)-6.0*  W85/lth/5.0 

Mas(5,l  l)=Mas(5,l  l)-60.0*W96/(lth**3)/5.0 

Mas(5, 1 2)=ALPHA*  W97-(W95-W86)/lth+W85/l  0.0+6.0*  W96/(lth*  *2) 

Mas(6,6)=(4.0*W66-W65*lth+2.0/15.0*W55*(lth**2))/lth 

Mas(6,7)=-Mas(5,8) 

Mas(6,8)=(60.0*W66-W55*(lth**2))/30.0/lth 

Mas(6,9)=ALPHA*(-W97)-(W95-W86)/lth+W85/l  0.0+6.0*  W96/(lth*  *2) 

Mas(6,12)=2.0*W96/lth-1.0/30.0*W85*lth+(W86-W95)/2.0 

Mas(6, 1 1  )=(W95-W86)/lth-ALPHA*  W97-W85/1 0.0-6.0*  W96/(Ith**2) 

Mas(6,10)=4.0*W96/lth+4.0/30.0*W85*lth-(W86+W95)/2.0 

Mas(7,7)=(ALPHA**2)*W77*lth+W55*6.0/5.0/lth+W66*12.0/(lth**3) 

Mas(7,7)=Mas(7,7)+2.0*W75*ALPHA 

Mas(7,8)=1.0/10.0*(10.0*W76*ALPHA-W55-W66*60.0/(lth**2)) 

Mas(7, 1 0)=W86/lth-ALPHA*  W97-W95/lth-W85/l  0.0-6.0*  W96/(lth*  *2) 

Mas(7,9)=ALPHA*(ALPHA*lth*W07+W05-W87)-6.0*W85/lth/5.0 

Mas(7,9)=Mas(7,9)-60.0*W96/(lth**3)/5.0 

Mas(7,ll)=ALPHA*(ALPHA*lth*W07+W05+W87)+6.0*W85/lth/5.0 

Mas(7,l  l)=Mas(7,l  l)+60.0*W96/(lth**3)/5.0 

Mas(7,12)=ALPHA*W97+(W95-W86)/lth-W85/10.0-6.0*W96/(lth**2) 

Mas(8,8)=(4.0*W66+W65*lth+2.0/15.0*W55*(lth**2))/lth 

Mas(8,9)=ALPHA*W97+(W95-W86)/lth+W85/10.0+6.0*W96/(Ith**2) 

Mas(8,10)=2.0*W96/lth-1.0/30.0*W85*lth-(W86-W95)/2.0 

Mas(8, 1 1  )=(W86-W95)/lth+ALPHA*  W97-W85/1 0.0-6.0*  W96/(lth**2) 

Mas(8,12)=4.0*W96/lth+4.0/30.0*W85*lth+(W86+W95)/2.0 

Mas(9,9)=lth*(ALPHA**2)*W00-2.0*ALPHA*W08+12.0*W99/(lth**3) 

Mas(9,9)=Mas(9,9)+6.0/5.0*W88/lth 

Mas(9, 1 0)=W99*6.0/(lth*  *2)-W09*  ALPHA+W88/1 0.0 

Mas(9,ll)=W00*(ALPHA**2)*lth-W99*12.0/(lth**3)-W88*6.0/5.0/lth 

Mas(9,ll)=Mas(9,ll) 

Mas(9,12)=W09*ALPHA+W88/10.0+6.0*W99/(lth**2) 

Mas(l  0,1 0)=4.0*  W99/lth-W98+2.0/l  5.0*  W88*lth 
Mas(10,ll)=-Mas(9,12) 

Mas(l  0, 1 2)=2.0*  W99/lth-W88*lth/30.0 

Mas(l  1,1  l)=ALPHA*(2.0*W08+ALPHA*lth*W00)+12.0*W99/(lth**3) 
Mas(l  1,1  l)=Mas(l  1,1  l)+6.0/5.0*W88/lth 
Mas(ll,12)=W09*ALPHA-6.0*W99/(lth**2)-W88/10.0 
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Mas(12,12)=W98+4.0*W99/lth+2.0/15.0*W88*lth 

Mass(l,l)=Mas(l,l) 

Mass(l,7)=Mas(l,2) 

Mass(  1 ,4)=Mas(  1 ,3) 

Mass(  1 , 1 0)=Mas(  1 ,4) 

Mass(l,2)=Mas(l,5) 

Mass(l,3)=Mas(l,6) 

Mass(  1 ,8)=Mas(  1 ,7) 

Mass(l  ,9)=Mas(l  ,8) 

Mass(l,5)=Mas(l,9) 

Mass(l,6)=Mas(l,10) 

Mass(l,ll)=Mas(l,ll) 

Mass(l,12)=Mas(l,12) 

Mass(7,7)=Mas(2,2) 

Mass(4,7)=Mas(2,3) 

Mass(7,10)=Mas(2,4) 

Mass(2J)=Mas(2,5) 

Mass(3,7)=Mas(2,6) 

Mass(7,8)=Mas(2,7) 

Mass(7,9)=Mas(2,8) 

Mass(5,7)=Mas(2,9) 

Mass(6,7)=Mas(2, 1 0) 

Mass(7,ll)=Mas(2,ll) 

Mass(7,12)=Mas(2,12) 

Mass(4,4)=Mas(3,3) 

Mass(4, 1 0)=Mas(3 ,4) 

Mass(2,4)=Mas(3,5) 

Mass(3,4)=Mas(3,6) 

Mass(4,8)=Mas(3,7) 

Mass(4,9)=Mas(3,8) 

Mass(4,5)=Mas(3,9) 

Mass(4,6)=Mas(3,10) 

Mass(4,ll)=Mas(3,ll) 

Mass(4,12)=Mas(3,12) 

Mass(10,10)=Mas(4,4) 

Mass(2, 1 0)=Mas(4,5) 

Mass(3 , 1 0)=Mas(4,6) 

Mass(8,10)=Mas(4,7) 

Mass(9,10)=Mas(4,8) 

Mass(5,10)=Mas(4,9) 

Mass(6,10)=Mas(4,10) 

Mass(  10,11  )=Mas(4, 11) 

Mass(10,12)=Mas(4,12) 

Mass(2,2)=Mas(5,5) 

Mass(2,3)=Mas(5,6) 


Mass(2,8)=Mas(5,7) 

Mass(2,9)=Mas(5,8) 

Mass(2,5)=Mas(5,9) 

Mass(2,6)=Mas(5 , 1 0) 

Mass(2,ll)=Mas(5,ll) 

Mass(2,12)=Mas(5,12) 

Mass(3,3)=Mas(6,6) 

Mass(3,8)=Mas(6,7) 

Mass(3,9)=Mas(6,8) 

Mass(3,5)=Mas(6,9) 

Mass(3,6)=Mas(6,10) 

Mass(3,ll)=Mas(6,ll) 

Mass(3,12)=Mas(6,12) 

Mass(8,8)=Mas(7,7) 

Mass(8,9)=Mas(7,8) 

Mass(5,8)=Mas(7,9) 

Mass(6,8)=Mas(7,10) 

Mass(8,ll)=Mas(7,ll) 

Mass(8,12)=Mas(7,12) 

Mass(9,9)=Mas(8 , 8) 

Mass(5,9)=Mas(8,9) 

Mass(6,9)=Mas(8,10) 

Mass(9,l  l)=Mas(8,l  1) 

Mass(9,12)=Mas(8,12) 

Mass(5,5)=Mas(9,9) 

Mass(5,6)=Mas(9, 1 0) 

Mass(5,ll)=Mas(9,ll) 

Mass(5,12)=Mas(9,12) 

Mass(6,6)=Mas(  10,10) 

Mass(6,ll)=Mas(10,ll) 

Mass(6,12)=Mas(10,12) 

Mass(l  1,1  l)=Mas(l  1,11) 

Mass(ll,12)=Mas(ll,12) 

Mass(12,12)=Mas(12,12) 

D0  51,R0W=2,12 
IM1=R0W-1 
DO  51,  COL=  1,IM1 
Mass(row,col)=Mass(col,row) 

5 1  continue 
end 

SUBROUTINE  Mass3  (rl,r2,hl,h2,lth,Mass) 
INTEGER  ROW,  COL 
REAL  Mass(12,12),Ith 

Al=(rl*hl/10.0)*(26.0*(lth**2)  +  7.0*(hl**2))/lth/7.0 


A2  =  (r2*h2/10.0)*(26.0*(lth**2)+7.0*(h2**2))/lth/7.0 
B1  =(rl*hl/10.0)  *  (11.0/21.0*(lth**2)+(hl**2)/12.0) 
B2  =  (r2*h2/10.0)  *  (11.0/21.0*(lth**2)+(h2**2)/12.0) 
Cl  =  (rl*hl/10.0)*(-9.0*(lth**2)+7.0*(hl**2))/lth/(-7.0) 
C2  =  (r2*h2/10.0)*(-9.0*(lth**2)+7.0*(h2**2))/lth/(-7.0) 
D1  =  (rl*hl/10.0)*(-13.0/42.0*(lth**2)  +  (hl**2)/12.0) 
D2  =  (r2*h2/10.0)  *  (-13.0/42.0*(lth**2)  +  (h2**2)/12.0) 
FI  =  (rl*hl/10.0)  *  (2.0/21. 0*(lth**3)+(hl**2)*lth/9.0) 
F2  =  (r2*h2/10.0)  *  (2.0/21. 0*(lth**3)+(h2**2)*lth/9.0) 
Z1  =  (rl*hl/10.0)  *  ((lth**3)/(-14.0)-lth*(hl**2)/36.0) 
Z2  =  (r2*h2/10.0)  *  ((lth**3)/(-14.0)-lth*(h2**2)/36.0) 
U1  =(rl*hl*lth/3.0) 

U2  =  (r2*h2*lth/3.0) 

Mass(l,l)=Ul 
Mass(l,7)=Ul/2.0 
Mass(7,7)=Ul 
Mass(4,4)=U2 
Mass(4,10)=U2/2.0 
Mass(10,10)=U2 
Mass(2,2)=Al 
Mass(2,3)=Bl 
Mass(2,8)=Cl 
Mass(2,9)=Dl 
Mass(3,3)=Fl 
Mass(3,8)=-Dl 
Mass(3,9)=Zl 
Mass(8,8)=Al 
Mass(8,9)=-Bl 
Mass(9,9)=Fl 
Mass(5,5)=A2 
Mass(5,6)=B2 
Mass(5,ll)=C2 
Mass(5,12)=D2 
Mass(6,6)=F2 
Mass(6,l  1)=-D2 
Mass(6,12)=Z2 
Mass(ll,ll)=A2 
Mass(ll,12)=-B2 
Mass(12,12)=F2 
DO  52,  ROW=2,12 
IMl=ROW-l 
DO  52,  COL=  1,IM1 
Mass(row,col)=Mass(col,row) 

52  continue 
end 


APPENDIX  3 

Ansys  Program 


/PREP7 

/TITLE,  SANDWICH  BEAM  FREE  VIBRATION,  Non-damping  study 

KAN, 2 

ET,I,82 

MP,NUXY,I,.3 

MP,EX,  1,29000000 

MP,DENS,1,. 000734 

N,  1,0,0 

N,2,0,.125 

N,3,0,.25 

N,4,0,.3125 

N,5,0,.375 

N,6,0,.4375 

N,7,0,.5 

N,8,0,.625 

N,9,0,.75 

N,10,.238486842,0 

N,ll,.238486842,.25 

N,12,.238486842,.375 

N,13,.238486842,.5 

N,14,.238486842,.75 

NGEN,20, 14,1,1 4„.476973684 

E,l,15, 17,3,10,16,1 1,2 

ENGEN,4,19,14,1 

EN,4, 7,21,23,9,13,22,14,8 

ENGEN,4,19,14,4 

ET,2,82 

MP,EX,2,345.8217 

MP,DENS,2,.000115 

MP,NUXY,2,.49 

MAT,2 

EN,2,3, 17,19,5,11, 18,12,4 
EN,3,5, 19,21,7,12,20, 13,6 
ENGEN,4,19,14,2 
ENGEN,4,19,14,3 
D,1,UY 

D,267,UX,„275 

FINISH 
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/SOLU 

MODOPT,SUBSP,5 

SOLVE 

FINISH 

/POSTl 

SET, LIST 

FINISH 


